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Abstract 



Gravitational and electrostatic interactions are fundamental examples of sys- 
tems with long-range interactions. Equilibrium properties of simple models 
with long-range interactions are well understood and exhibit exotic behaviors : 
negative specific heat and inequivalence of statistical ensembles for instance. 

The understanding of the dynamical evolution in the case of long-range in- 
teracting systems still represents a theoretical challenge. Phenomena such as 
out- of- equilibrium phase transitions or quasi- stationary states have been found 
even in simple models. 

The purpose of the present thesis is to investigate the dynamical properties of 
systems with long-range interactions, specializing on one-dimensional models. 
The appropriate kinetic description for these systems is the Vlasov equation. 
A statistical theory devised by D. Lynden-Bell is adequate to predict in some 
situations the outcome of the dynamics. A complementary numerical simulation 
tool for the Vlasov equation is developed. 

A detailed study of the out-of-equilibrium phase transition occuring in the 
Free-Electron Laser is performed and the transition is analyzed with the help 
of Lynden-Bell's theory. Then, the presence of stretching and folding in phase 
space for the Hamiltonian Mean-Field model is studied and quantified from 
the point of view of fluid dynamics. Finally, a system of uncoupled pendula 
for which the asymptotic states are similar to the ones of the Hamiltonian 
Mean-Field model is introduced. Its asymptotic evolution is predicted via both 
Lynden-Bell's theory and an exact computation. This system displays a fast 
initial evolution similar to the violent relaxation found for interacting systems. 
Moreover, an out-of-cquilibrium phase transition is found if one imposes a self- 
consistent condition on the system. 

In summary, the present thesis discusses original results related to the oc- 
curence of quasi-stationary states and out-of-equilibrium phase transitions in ID 
models with long-range interaction. The findings regarding the Free-Electron 
Laser are of importance in the perspective of experimental realizations of the 
aforementioned phenomena. 
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Resume 



Les interactions gravitationnelles et electrostatiques sont deux exemples fonda- 
mentaux de systemes en interaction de longue portee. Les proprietes d'equilibre 
de modeles simples en interaction de longue portee sont bien comprises et 
revelent des comportemens exotiques : capacite specifique negative et inequivalence 
des ensembles statistiques par exemple. 

La comprehension dc revolution dynamiquc dans le cas de systemes en inter- 
action de longue portee represente encore actuellement un defi theorique. Des 
modeles simples presentcnt des proprietes telles que des transitions de phase 
hors d'equilibre ou des etats quasi- stationnaires. 

Le but de la presente these est d'etudier les proprietes dynamiques de systemes 
en interaction de longue portee pour des modeles a une dimension. La de- 
scription cinctique adequate est donnee par l'equation de Vlasov. Une theorie 
statistique proposee par D. Lynden-Bell est appropriee pour predire dans cer- 
taines situations l'aboutissement de la dynamique. Un outil de simulation pour 
l'equation de Vlasov complete cette approche. 

Une etude detaillee de la transition de phase dans le Laser a Electrons Libres 
est presentee et la transition est analysee a l'aide de la theorie de Lynden-Bell. 
Ensuite, la presence d'etirement et de repliement est etudiee dans le modele 
Hamiltonian Mean-Field cn analogie avec la dynamique des fluides. Enhn, un 
systeme de pendules decouples dont les etats asymptotiques sont similaires a 
ceux du modclc Hamiltonian Mean-Field est introduit. Son evolution asympto- 
tique est predite par la theorie de Lynden-Bell et par une approche exacte. Ce 
systeme presente une evolution initiale rapide similaire a la relaxation violente 
presente dans des modeles plus compliques. De plus, une transition de phase 
hors d'equilibre est trouvee si une condition d'auto-consistence est imposee. 

En resume, la presente these comporte des rcsultats originaux lies a la 
presence d'etats quasi-stationnaires et de transitions de phase hors d'equilibre 
dans des modeles unidimcnsionnels en interaction de longue portee. Les resultats 
concernant le Laser a Electrons Libres offrent une perspective de realisation 
experimentale des phenomenes decrits dans cette these. 
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Chapter 1 

Introduction 



The long-range character of the interaction modifies significantly the macroscopic 
behaviour of physical systems. We present general properties of such systems. The 
reasons motivating their investigations are observations in the context of galactic dy- 
namics and in systems of charged particles. We specialize on the particular topic of 
the collisionless evolution in ID models of systems with long-range interactions. The 
Vlasovian description of these systems has been the focus of recent interest, inducing 
the need for a numerical implementation. 



1.1 Motivations 

Gravitational and electrostatic interactions are fundamental examples of sys- 
tems with long-range interactions. The long-range character of the interaction 
has a profound impact on the thermodynamical and dynamical features of these 
systems. A treatment for neutral plasmas allows one to introduce a screening 
while for gravitational systems, non-neutral plasmas and beams of charged par- 
ticles the unshielded interactions manifest long-range properties. 

Despite its commonness in physics, the detailed study of long-range interact- 
ing systems has only recently attracted important interest from the statistical 
physics community, and has led to the organization of conferences and sum- 
mer schools (Dauxois et al., 2002; Campa et al., 2008; Dauxois et al., 2008). 
The focus of the long-range community has been mainly devoted to toy models 
and equilibrium statistical mechanics in the microcanonical ensemble, bring- 
ing out peculiarities such as ensemble inequivalence, negative specific heat and 
non-additivity. These phenomena do not exist in systems with short-range in- 
teractions which, in contrast, are additive. 

The aforementioned meetings brought together specialists from different 
fields, experimentalists and theoreticians alike, with the objective of bringing 
links between theory and experiments in the coming years. This perspective 
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seems especially promising for trapped plasmas (Drewsen ct al., 2008), dipolar 
interactions (Menotti et al., 2008) or Free-Electron Lasers (Cur bis et al., 2007). 

Another trend is to focus on the dynamical phenomena found in systems with 
long-range interactions. More precisely, the slow relaxation to equilibrium, or 
the lack of relaxation, has attracted much interest. This feature offers challenges 
and motivations in the context of gravitational systems or wave-particle inter- 
actions. The transitory stage of the dynamics in which the system may remain 
trapped for long times before eventually reaching thermodynamical equilibrium 
is referred to as a Quasi-Stationary State (QSS) (Antoni and Ruffo, 1995; An- 
toniazzi et al., 2007b). The existence of out-of-equilibrium phase transitions 
separating different dynamical regimes contributes to enrich the interest in dy- 
namical studies (Antoniazzi ct al., 2007c; Antoniazzi et al., 2007b; Latora et al., 
1998; de Buyl et al., 2009). 

A precise description of the dynamical evolution can be achieved thanks to 
kinetic theory. Chavanis has dedicated a series of papers ((Chavanis, 2006a) to 
(Chavanis, 2008c)) to the application of kinetic theory to systems with long- 
range interactions. In addition, a number of articles have focused on the Vlasov 
equation (Yamaguchi et al., 2004; Barre et al., 2004; Barre et al., 2006; Antoni- 
azzi et al., 2007a), which is adequate to describe systems where collisional effects 
are absent. This equation is also called collisionless Boltzmann equation in as- 
trophysics. The Vlasov equation allowed the understanding of stability and 
dynamical properties, and is the basis for Lynden-Bell's theory, a statistical 
theory aimed at the prediction of collisionless evolution following the so-called 
violent relaxation. The theory was introduced in Ref. (Lynden-Bell, 1967). 

The purpose of this thesis is to study the dynamical properties of the Vlasov 
equation, especially in the process of relaxation. We consider the Hamiltonian 
Mean-Field model (Antoni and Ruffo, 1995) which can be viewed as a paradigm 
to study systems with long-range interactions and the Colson-Bonifacio model 
for the Free-Electron Laser (Bonifacio et al., 1986), an example of wave-particle 
interaction. We study the Vlasov equation itself through direct numerical res- 
olution, going beyond molecular dynamics simulations. We detail a numerical 
procedure to perform these computations. We finally use an analogy with a set 
of uncoupled pendula to refine our understanding of the dynamics. 

1.2 Examples of systems with long-range inter- 
actions 

We give in this section a number of physical situations where the long-range 
nature of the interaction needs to be taken into account. A number of them 
derives from the fact that charged particles interact via a 1/r potential 1 , and 
in some situations, only the mathematical description is effectively long ranged 
while the underlying physical system is not. 



In this introduction, r is generically considered the distance between particles. 
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Charged particles. Given a density of particles p, the interaction potential 
$ between particles under Coulombian interaction is given by the solution of 
Poisson's equation : 

A$ = -p (1.1) 

where A is the Laplacian operator. The solution is an inter-particle potential 
V(r) of 1/r in 3D, — logr in 2D and — r in ID. These potentials are considered 
long range : the decay of V(r) in function of r is weak. They correspond to 
point particles in 3D, infinite charged parallel wires in 2D and sheets in ID. The 
infinite charged wires are an adequate description for beams of charged particles, 
which we describe in a separate paragraph. The ID system is equivalent to a set 
of parallel charged sheets (see Fig. 1.1) and does not correspond to a physically 
observable system. It has served as a theoretical model to study plasma physics 
since a long time [see for instance Ref. (Dawson, 1983)]. 




Figure 1.1: A system of parallel sheets. The position is given along the direction 
normal to the sheets. 



Models of gravitational interaction. The relaxation time of galaxies, i.e. 
the time after which we expect the galaxy to have reached thermodynamical 
equilibrium, can be as long as the age of the universe ( Chandrasekhar, 1942). 
Galactic dynamics thus needs to be analyzed from a collisionless perspective. 
The ID sheet model presented belows allows for a simplified study of the phe- 
nomenon. In this thesis, the topic of gravitational interaction motivates the use 
of the Hamiltonian Mean-Field model defined in chapter 3 : it can be obtained 
by retaining only the first Fourier mode of the ID sheet model. 

Given a density of particles p, the gravitational potential $ between particles 
is given by the solution of Poisson's equation : 

A0 = P (1.2) 
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where A is the Laplacian operator. The solution is an inter-particle potential 
V(r) of — 1/r in 3D while in ID the solution is r. 

The ID sheet model has been the subject of numerous numerical studies 
(Hohl, 1968; Luwel et al., 1984; Luwel and Severne, 1985) motivated by its 
simplicity and the existence of exact (up to machine precision) algorithms to 
simulate its dynamics. It is usually called the gravitating sheets model, because 
the ID version corresponds to a hypothetical system of parallel rigid sheets, 
described by their positions along the normal axis (see Fig 1.1). 

Beams of charged particles. In a beam of charged particles (ions or elec- 
trons), the longitudinal velocity along the beam propagation direction (say z) 
can be considered a constant. The particles can then be described in the x — y 
plane only, with a logarithmic repulsive interaction. 

This system has been the subject of experimental work on vortex dynamics 
and turbulence (Driscoll and Fine, 1990), and is also a relevant description for 
charged particles circulating in an accelerator (Benedetti et al., 2006). 

Wave-particle interactions. Wave-particle interactions are useful to de- 
scribe situations in which the inter-particle forces can be neglected in compar- 
ison with the mean-field coupling provided by an electromagnetic wave. While 
in this thesis we interpret this model as a description of a Free-Electron Laser 
device (see chapter 5 for details), it is of more general interest in plasma physics 
(Elskens and Escande, 2003; del Castillo-Negrete, 2002). 

The coupling via a wave allows for effective long-range interactions and pro- 
vides rich dynamical properties. The reduction of the wave to one component al- 
lows for a simplified treatment. The Colson-Bonifacio model that we investigate 
in this thesis is however relevant to describe features of realistic Free-Electron 
Laser devices (Bonifacio et al., 1990; Bonifacio et al., 1986; Curbis et al., 2007). 

Other systems of interest. Two-dimensional flows can be considered effec- 
tively long-ranged. The vorticity field — [where u x (u y ) is the velocity in 
the x (y) direction] displays a logarithmic interaction. A statistical theory sim- 
ilar to Lynden-Bell's theory predicts the evolution of the flow (Chavanis et al., 
1996). 

The dipolar interaction is also long-range, going as 1/r 3 in 3D. The dipolar 
interaction is a candidate for experiments on long-range phenomena (Menotti 
et al., 2008). 

We call "Small system" a system in which the range of the interaction is of 
the order of the size of the system. In that setup, each particle feels the presence 
of the complete set of particles in the system. Thermodynamics may be affected 
by this property and some features such as non-additivity or negative specific 
heat may appear (Chomaz and Gulminclli, 2002). 
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1.3 Thermodynamical behaviour of long-range 
interacting systems 



The range of the interaction potential defines a neighborhood for each particle. 
In order to model the properties of an ensemble of particles, we take into account 
the coupling of a particle to all of its neighbours. As an example to illustrate 
these properties, we take the case of the Ising model, the Hamiltonian of which 
is the following: 

N 

H = -J J2 S i S 3 ( L3 ) 
<i,j> = l 

where J is a coupling constant, Si = —1,1 the spin value at site i and the 
brackets in the sum means that the sum is on nearest neighbours only. In the 
simple setup of a square lattice, a spin has 2 nearest neighbors in one dimension, 
4 in two dimensions and 6 in three dimensions. It is clear that the interaction 
energy for a single site is bounded by its number of neighbours. 

Extending the range of the interaction to infinity leads to mean-field models: 

N / N \ 2 j 

H = -j'£S i S J = -U^S i ) =- J -NW (1.4) 

i<j=l \i=l / 

where the sum runs over all pairs i,j and 

M=l]>> (1.5) 

i 

is the magnetization. The bound of the energy at each site is now growing as 
N and the global energy scales as iV 2 . While H can be made extensive by an 
appropriate rescaling, it cannot be made additive. 



Lack of additivity. In systems with short-range interactions, the interactions 
between two subsystems take place at the interface between the two subsystems. 
The bulk energy then grows with the volume, while the interfacial interaction 
only grows in proportion to the surface. The interfacial term thus become 
negligible in the thermodynamic limit. 

In long-range interacting systems, the interfacial interaction needs to take 
into account all particles and therefore grows in proportion to the volume. It 
cannot be neglected anymore when compared to the bulk energy. We illustrate 
that property using systems (1.3) and (1.4). 

Let us first consider the Hamiltonian (1.3) in two dimensions. If we divide 
the systems into two parts a and b whose spins are labeled respectively Sf and 
Sf. When the subsystems are placed next to each other (see Fig. 1.2), the total 
energy H a+b can be expressed as: 

H a+b = H a + H b_ jJ2S?S b j (1.6) 
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where the sum runs over nearest neighbours only (here along the boundary 
indicated in Fig. 1.2). The surface term scales as vN and therefore becomes 
negligible compared to H a+b oc N in the thermodynamic limit. Namely, we 
may write H a+b w H a + H b . 



+ + + + + + 
+ + + + + + 
+ + + + + + 

:+"+ + + + + 



Figure 1.2: Additivity of a system 
of spins interacting with their near- 
est neighbours : The division of the 
system into two parts is such that 
the overall energy is approximately 
the sum of the energy of the two 
parts considered separately. The dot- 
ted line indicates the region that 
contributes to the difference between 
H a + H b and H a+b ', negligible with 
respect to the size of the bulk contri- 
bution. 



We now turn to the system (1.4). Dividing it into two subsystems, we write 
the energy of the total system as: 

N a N b N a N b 

H a +b = S a + J- S b f = H a + H b - JC£ Sf)(£ S b ) (1.7) 

i— 1 i—1 i—1 i—1 

where the non-additive part cannot be neglected anymore (see Fig. 1.3). 



Ensemble inequivalence. The usual equivalence between the microcanon- 
ical and canonical ensembles can be proved on general grounds for systems 
satisfying the additivity property. In systems with long-range interactions, the 
equivalence of ensembles is no longer granted and needs to be checked individ- 
ually for model under consideration. In the case of inequivalence, the adequate 
physical description depends on the constraints put on the system: isolated sys- 
tems correspond to the microcanonical ensemble while systems in contact with 
a heat bath correspond to the canonical ensemble. 

Ensemble equivalence holds for the models considered in this thesis. We refer 
the reader to Refs. (Barre et al., 2001; Barre, 2003; Barre et al., 2005; Campa 
et al., 2009) to find a discussion of ensemble inequivalence. Reference (Barre 
et al., 2001) investigates the phase diagram of the Blume-Emery-Grifhths model 
with infinite range interactions in which ensemble inequivalence is found. 
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Figure 1.3: Non-additivity of a mean- 
field spin system. The division of the 
system into two parts cannot approxi- 
mate the properties of the full system. 
The dotted line indicates the region 
that contributes to the difference be- 
tween H a + H b and H a+b . In the case 
of long-range interaction, that region 
can be of the size of the complete sys- 
tem. 



Negative specific heat. The positive character of the specific heat also fol- 
lows from additivity. A negative specific heat can be found in the microcanonical 
study of systems with long-range interactions. The occurrence of negative spe- 
cific heat implies the property of ensemble inequivalence : it can be proved to 
be positive in the canonical ensemble (Barre et al., 2001; de Buyl et al., 2005). 



1.4 Dynamical properties 

Numerical experiments on the Hamiltonian Mean-Field model (HMF) have led 
to the discovery of a particular dynamical phenomenon, the so-called quasi- 
stationary states (QSS). Starting from an initial condition, the system evolves 
in time towards a QSS before eventually reaching equilibrium. For a time which 
turns out to diverge when the number of particles is increased, the macroscopic 
quantities characterizing the system remain distinct from the value expected at 
equilibrium. Situations of interest arise when the time within which we observe 
the system is such that equilibrium statistical mechanics does not provide a 
good prediction. It is consequently a challenge to understand these QSS and 
ultimately to predict their properties. 

An unusual feature is that the emergence of dynamical properties in systems 
with long-range interactions greatly depends on the initial condition, while the 
initial condition is usually considered irrelevant in classical statistical mechan- 
ics. It eventually leads to an out- of- equilibrium phase transition, a change in 
the macroscopic observable found when the system is kept away from thermo- 
dynamical equilibrium. 
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Quasi-Stationary States (QSS). Quasi-Stationary States (QSS) arc dy- 
namical states in which the system displays features of equilibrium but that 
does not correspond to the equilibrium predicted by statistical mechanics in the 
following sense: 

• The system eventually experiences a change in its macroscopic variable. 

• That change takes place at a time that increases when the number of 
particles considered increases. 

• The macroscopic quantities characterizing the system do not correspond 
to those predicted by statistical mechanics. 

Their observation was reported in numerical experiments and have since at- 
tracted a lot of attention in the literature (Yamaguchi et al., 2004; Barre et al., 
2006; Chavanis, 2006d; Antoniazzi et al., 2007b). The time spent in a QSS may 
be very large, and a prediction of that state is highly desirable. 

The Vlasov equation allows one to cast a theoretical framework to under- 
stand the QSSs. The time of validity of the description by the Vlasov equation 
is investigated in (Jain et al., 2007) and discussed in relation with numerical 
simulations in (Antoniazzi et al., 2007a). 

A statistical theory devised by Lynden-Bell and appropriate for the Vlasov 
equation (Lynden-Bell, 1967) has been used with success to predict the QSS for 
the Hamiltonian Mean-Field (Antoniazzi et al., 2007b; Antoniazzi et al., 2007c), 
and to describe them in the Free-Electron Laser (Barre et al., 2004; Curbis et al., 
2007; de Buyl et al., 2009) 

Numerical simulations on gravitational systems report cases of slow relax- 
ation (Luwel et al., 1984). Yamaguchi applied a slightly modified Lynden-Bell 
theory to the ID sheet model and discussed the concept of QSS in that context 
(Yamaguchi, 2008). 

Core-Halo distribution. A core-halo distribution is a non-monotonous par- 
ticle distribution. It has been first observed in the Gravitating Sheets Model 
(Yamashiro et al., 1992), but can be found also in the Hamiltonian Mean-Field 
model. No statistical theory has been predicting these distributions which are 
understood as a purely dynamical phenomenon. 

Out-of-equilibrium phase transition. An out-of-equilibrium phase transi- 
tion is a dynamical phenomenon that can cause a system to evolve into radically 
different regimes, depending on the parameters of the initial condition. In the 
Hamiltonian Mean-Field model, for instance, such a transition is found numer- 
ically and can be predicted by Lynden-Bell's theory (Antoniazzi ct al., 2007c). 
Further analysis revealed that the transition is related to the occurrence of reg- 
ular structures in phase space, from a one-cluster state to a two-clusters state 
(Bachclard et al., 2008). In the Free-Electron Laser, a similar out-of-equilibrium 
phase transition exists, and the analysis of the associated phase space structures 
is expected to display a rich phenomenology (de Buyl et al., 2009). 
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The Vlasovian approach. The description of the dynamical evolution by 
the Vlasov equation is valid in the thermodynamic limit A — > oo. Its use 
helps the understanding of dynamical phenomena for systems with long-range 
interactions. The description of the quasi-stationary states of the dynamics by 
Lynden-Bell's theory, which is based on the Vlasov equation, is a remarkable 
theoretical achievement. It has been to predict successfully the magnetization in 
the Hamiltonian Mean-Field model (Antoniazzi et al., 2007c; Antoniazzi et al., 
2007b) and the intensity in the free-electron laser (Barre et al., 2004; Curbis 
et al., 2007). 

The kinetic description the Vlasov equation provides does not take into ac- 
count the number of particles A. However, the numerical simulations performed 
in the literature are almost exclusively of the A-body type. An exception is 
given by Ref. (Antoniazzi et al., 2007a). The impact of the number of particles 
on the dynamics has been quantified, with special interest on the time scale for 
which the Vlasov equation and the A-body description coincide (Yamaguchi 
et al., 2004; Jain et al., 2007). 

This thesis investigates the dynamical evolution of the Vlasov equation 
through direct numerical simulation in which the number of particles A is ir- 
relevant. This approach brings a new tool to the study of models with long- 
range interactions and helps to shed light on the detailed phase space dynamics 
(dc Buyl, 2009). 

Relaxation scenario The authors of Ref. (Barre et al., 2006) propose that 
the relaxation scenario in systems with long-range interactions follows the fol- 
lowing steps : a first stage is called violent relaxation, its time scale does not 
depend on the number of particles A ; following violent relaxation, the system 
may get trapped in a quasi-stationary state (QSS) and eventually attains equi- 
librium after a time that depends on the number of particles. This scenario is 
depicted in Fig. 1.4 The dependence on the number of particles present in the 
system is studied in Refs. (Yamaguchi et al., 2004) and (Jain et al., 2007). 

The use of direct simulations of the Vlasov equation allows one to work in the 
limit A — > oo and provides an useful addition to A-body simulations regarding 
the lifetime of the QSSs. Vlasov simulations also offer the possibility to study 
the stability of arbitrary solution of Lynden-Bell's theory (de Buyl et al., 2009). 
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Figure 1.4: The relaxation scenario for Vlasov dynamics, as proposed in 
Ref. (Barre et al., 2006). An initial condition evolves through violent relax- 
ation on a timescale of order 0(1). If the dynamics gets trapped into a quasi- 
stationary state (QSS), the time needed to reach equilibrium increases with the 
number of particles : t « N s or t w log N (Jain et al., 2007). 
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1.5 Outline of the thesis 



Chapter 2 is a general introduction to the Vlasov equation and its properties. 

Chapter 3 gives recalls briefly the definition of the Hamiltonian Mean-Field 
model and presents some of its general properties. 

Chapter 4 details the numerical procedure required to solve the Vlasov equa- 
tion. A study of the Hamiltonian Mean-Field model is performed and the limi- 
tations of the method are discussed. The case of the Free-Electron Laser is also 
given. 

Chapter 5 introduces a model for the Free-Electron Laser, its associated 
Vlasov equation and presents the study of the out-of-cquilibrium phase transi- 
tion. 

Chapter 6 studies the relaxation in the Hamiltonian Mean-Field model on 
the basis of numerical simulations of the Vlasov equation. The analysis benefits 
of an analogy with fluid dynamics. 

Chapter 7 concerns a set of uncoupled pendula. In that model, we derive 
an exact solution to the asymptotic evolution of the Vlasov dynamics under an 
ergodic hypothesis. 

Chapter 8 concludes this thesis and offers some perspectives to its results. 

The table of contents contains a short abstract of each chapter and consti- 
tutes a good starting point to the reading of this thesis. 



27 



28 



Chapter 2 

The Vlasov equation 



We introduce the Vlasov equation as the continuum limit in systems with long- 
range interactions. We discuss its derivation via the BBGKY hierarchy of kinetic 
theory and the conservation properties of Vlasov dynamics. We present Lynden-BelPs 
theory for the statistical mechanics of violent relaxation. 



The dynamical evolution of particles under mutual interactions is described, 
within the range of validity of classical mechanics, by Newton's equation, ruling 
the evolution of each particle's position and velocity. When the number of 
particles goes to infinity, it is impossible to follow the particles' coordinates 
individually. Instead, we will be interested in the probability with which we 
find the set of particles at given positions and velocities. The knowledge of 
that probability is sufficient to compute the average of macroscopic quantities. 
The equation ruling the evolution of that probability is called the Liouville 
equation, and its resolution is as complicated as solving Newton's equations 
for the ensemble of individual trajectories, unless simplifying assumptions can 
be made. The objective of kinetic theory is to identify physically interesting 
approaches that allow such simplifications and to understand the statistical 
behaviour of dynamical phenomena (Balescu, 1997). For instance, we restrict 
our attention in this thesis to the probability of finding one particle at a given 
position and velocity, irrespective of the remaining particles. 

This chapter deals with one approach, namely the Vlasov equation, whose 
domain of validity is interesting for systems with long-range interactions which 
are considered in this thesis. This chapter is organized as follows: we discuss 
the phase space description in section 2.1, introduce the Vlasov equation in 
section 2.2 and finally present Lynden-Bell's theory in section 2.3 
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2.1 The phase space description 



We limit ourselves, as this is not a limitation for this thesis, to the study of 
periodic systems in which the position is labeled by an angle 6 and is defined in 
the interval [— 7r; n]. The associated momentum is denoted p. 

2 

Given a Hamiltonian H for particles with a kinetic energy and an inter- 
particle potential energy V(r) (where r is the distance between two particles), 
the evolution of and p obeys the following equations: 

H = Ei + ^E^^i) 

. dH 

Si = ^— = Pi 
dpi 

* - -^ = -E^W-^l) (2-1) 

To display the state described by the system (2.1), we project all particles' 
positions and momenta in /x-space, the space of single-particle coordinates and 
momenta; we will call it equivalently phase space in this thesis. Fig. 2.1 gives 
an example of 200 particles uniformly distributed and moving to the right (i.e. 
with positive velocity). 



0.4 

0.2 
&, 0.0 
-0.2 
-0.4 

-3-2-10 1 2 3 

9 

Figure 2.1: An example of distribution of 200 point particles in phase space. The 
particles are uniformly distributed and move with a positive average velocity of 
about 0.3 . 

The system can be described alternatively by a distribution with support is 
phase space: / (0,P) [here, (©, P) belongs to the iV-particle phase space : it 
represents the full set {{0i,Pi)} of positions and velocities]. f N (<d,P) depends 
on the time t, but this dependence is not explicitly written to simplify the 
notation. In the Vlasov regime (see paragraph 2.2), we reduce the problem 
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to the study of the one-particle distribution function : Nf(0,p) d0 dp is the 
number of particles having their position and momentum contained in a cell of 
volume dO dp centered at (9,p). The distribution function f(9,p) is normalized 
to 1: 

J d6 dp f(9,p) = 1 (2.2) 

We reproduce the example of Fig. 2.1 with more particles in Fig. 2.2. In 
the thermodynamic or kinetic limit, the distribution function f(9,p) approaches 
a smooth distribution. We illustrate in Figs. 2.2 and 2.3 the change from the 
particles' point of view to the kinetic one. 




n, o.o 

-0.2 
-0.4 

-3-2-10 1 2 3 




e e 

Figure 2.2: (left panel) An example of distribution of 1500 point particles in 
phase space. The particles are uniformly distributed and move with a positive 
velocity of about 0.3 . 

Figure 2.3: (right panel) A kinetic description of the situation depicted in Fig. 
2.2 with a distribution function. The color code on the right indicates the value 
of the distribution function /. 



2.2 Kinetic equations 
2.2.1 Formalism 

We introduce the formalism of the reduced distribution functions (following 
(Balescu, 1997)). Given the probability F(Q, P)dQdP that a system of A par- 
ticles finds itself with positions and velocities in the region of volume dQdP 
around (0,P), where = {6\, . . . ,0/v} and P — {pi, . . . ,Pn}, we define the 
s-particles reduced distribution functions: 

A' f 

f s ((e 1 , Pl ),...,(9 s , Ps )) = (JV _ a) , J d9 s+1 d Ps+1 ..d9 N d PN F(G,P) (2.3) 

F is implicitly dependent on the time t. We will not write down the explicit 
time dependence of F, and of the subsequent /'s and <?'s defined in this section, 
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unless needed for clarity. The evolution of F is ruled by the Liouville equation: 

^F(@,P)=CF(Q,P) (2.4) 

where C is the Liouvillian operator, defined through a Poisson bracket: 

\^[dH_dF_ _ dH_dF_\ 

Using the definition (2.3), and computing the evolution equations for the f s 's 
leads to the BBGKY hierarchy of equations, equivalent to the Liouville equation 
(Balescu, 1997). 

Under the assumption that no correlation exists between particles, we may 
write the uncorrelated one-particle distribution function as : 

S 

/rm.pi). • ■ • . (°s,p s )) « n hWhPi)- ( 2 -6) 

3 = 1 

Beyond that assumption, one may use the correlation functions g s to quan- 
tify the deviation from this ideal behaviour, g 2 being defined by the following 
relation: 

/ 2 ((01,Pl),(02,P2)) =fl(01,Pl)fl(02,P2)+g2{{9l,Pl),{02,P2)). (2-7) 

Similar relations define g s for s > 3 in the cluster representation. 

The evolution of /i, deduced from the Liouville equation, is given by : 

iMOuPi) +pi^m,Pi) = 

Jd9 2 d P2 ^-V(\9 1 -9 2 \)-^- i f 1 (e 1 ,p 1 )f 1 (9 2 ,p 2 ) 

+ Jd9 2 d P2 ^V{\0! 9 2 \) ■ - ^) ff 2(((?l,Pl), (<?2,P2)) 

(2.8) 

The meaning of the terms appearing in Eq. (2.8) is the following: 

• ADV = P\-^f{9\,p\) is an advection term. It causes the displacement 

of molecules (or phase space elements in the kinetic description) according 
to their velocities. When only this term is present in the equation, the 
situation is called free streaming; each particle moves only at its constant 
velocity. 

. VLA = / d9 2 d P2 - 9 2 \) ■ ^-f 1 (e 1 ,p 1 )f 1 (e 2 ,p 2 ) is the Vlasov 

term. Given the interaction potential V and the distribution of particles 
in space, this term causes the particles to accelerate or decelerate. 
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. COL = /d0 2 d P2 ^-V^ - 9 2 \) ■ - ^)g 2 ((e 1 ,p 1 ),(6 2 , P2 )) is the 

collisional term. It causes the particles to accelerate or decelerate, as the 
Vlasov term, but depends on g 2 , implying that information is needed from 
a superior term in the BBGKY hierarchy (absence of closure). 

An external potential can easily be added to Eq. (2.8). It takes the following 
form [in the right hand side of Eq. (2.8)]: 

dV cxt df , nn . 

where V cxt (8) is a potential that does not depend on /. The study of a kinetic 
equation containing only such a term is given in chapter 7 



2.2.2 Boltzmann's equation 

Assuming the following conditions: 

1. A low density of particles ; 

2. A short-range interaction ; 

3. Weak inhomogeneity with respect to the range of the interaction ; 

4. Molecular chaos : the famous Stosszahlansatz ; colliding particles are con- 
sidered uncorrelated ; 

we obtain the Boltzmann equation, the most famous kinetic equation. It can be 
derived on the basis of a gain and loss methodology. The variation of fi(6i,pi) 
is given, aside from the advection, by the number of colliding particles resulting 
in a position and velocity (0\,pi) minus the number of particles undergoing 
collisions that make them leave the position (6\,pi). 
The Boltzmann equation takes the form: 

f + P f = G-L (2.0) 

where G represents the gain term and L the loss term, and G— L is the collisional 
term. 

The right hand side of Eq. (2.10) is the driving force of the relaxation to 
equilibrium. It represents binary encounters in the system, events that cause two 
nearby particles to change their momenta, under the assumption of Boltzmann's 
equation. An important property of Boltzmann's equation is the existence of 
an H-theorem proving the increase of entropy in the course of time and leading 
the system to a Boltzmann-Gibbs distribution. 
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2.2.3 The Vlasov equation 

In systems with long-range interactions, the aforementioned Vlasov term (see 
the term VLA in section 2.2.1) becomes non-negligible. 

Starting from Eq. (2.8) and assuming that the interaction potential is weak 
V(r) — 0(A), we establish the following estimation : 

f(6,p) = 0(\°) ;. 92 = 0(A 1 ) (2.11) 

Keeping the terms of Eq. (2.8) up to order 0(X 1 ), we obtain the Vlasov equation: 

where we have again dropped the i subscript. In this equation, 

V[f}(9) = J dff dp' f(6',p')V(\6 - 6'\) (2.13) 

is the self-consistently computed potential. As a result of the dependence of V 
on /, the Vlasov equation acquires a nonlinear character. 

Braun & Hepp (Braun and Hepp, 1977) proved the weak convergence of 
the A^-body dynamics to the Vlasov equation under the assumption of weak 
coupling. 

2.2.4 Properties of the Vlasov equation 

The Vlasov equation possesses the form of an advection equation in phase space. 
As a result, the levels of f(9,p) arc conserved by the dynamics: the volume of 
phase space in which / takes a value comprised between r] and ij + di] is constant. 
This property can be made explicit if we write the equation in a conservative 
form : 

i'-i'+'I'-sl'- (2 - 14) 

where 4z is the total time derivative. A major consequence is that an infinite 
number of quantities, the Casimirs, are conserved. For any function s, g|C s [/] = 
0, where 

C s [f] = J d9 dp s(f(e,p)) (2.15) 

Additionally, the preservation of all levels r\ of / implies that the maximum and 
the minimum of / remain constant. 

We define the Li norms, that we use in chapter 4, 

Lilf] = Jd9dp (f(9,p)Y , (2.16) 

where i is a positive integer. 
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Consider the one-particle Hamiltonian denned by : 

2 

H 1 (9,p) = ^ + V[f](9) (2.17) 

Any distribution depending only on Hi(6,p), f(8,p) = f(H\(0,p)), is a station- 
ary solution of the Vlasov equation. In particular, for homogeneous systems 
with no external force, this implies that any distribution of the momenta alone 
f(0,p) — ip(p) is stationary. The stability of these solutions needs to be studied 
individually; stationary stable solutions are candidates for equilibrium in Vlasov 
dynamics while this is not the case for unstable stationary solutions. 



2.2.5 Filament at ion 

Even in the simplest cases, the dynamical evolution of the Vlasov equation 
displays filamentation. Filamentation is the thinning of the structure in phase 
space. We illustrate filamentation in the motion of free streaming, given by the 
Vlasov equation with no force term. 

We display in Fig. 2.4 a waterbag initial condition 1 . A waterbag is a distri- 
bution that is constant inside a region of phase space, it is defined as : 

JMA- P ifbl<Ap, 
f(6,p)={ |0|<A0, (2.19) 

otherwise. 

The evolution of an initial / given by Eq. 2.19 is portrayed in Figs. 2.5 and 
2.6. Figure 2.5 shows the free-streaming initial evolution. Then, in Fig. 2.6 , 
we observe filaments. The width of the filaments and the space between them 
decreases in the course of time, and eventually becomes too small to allow a 
numerical description. 



2.3 The statistical theory of Lynden-Bell 

A feature of some systems with long-range interactions is that the time spent 
without reaching thermodynamical equilibrium may be large. It has been re- 
ported to diverge when the number of particles increases in numerical observa- 
tions (Barre et al., 2004; Yamaguchi et al., 2004). This result causes problems if 
one considers the thermodynamic limit N — > oo that is relevant for the Vlasov 
equation. For the time scale considered (experimental or observational), equi- 
librium statistical mechanics provides no useful prediction. In the framework of 

lr The waterbag initial condition is used here because of its frequent use in the following of 
this thesis. 
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Figure 2.4: (left): A waterbag initial condition. 

Figure 2.5: (right): Early evolution of the initial condition given in Fig. 2.4 by 
Eq. 2.18 Each horizontal level moves at a constant speed. 




Figure 2.6: Later evolution of the 
free-streaming. 







Vlasov dynamics, Lynden-Bell proposed (Lynden-Bell, 1967) a statistical theory 
to predict the result of the so-called violent relaxation process. 

Lyndon-Bell's theory, and more generally the use of the Vlasov equation, 
allows to understand the properties of iV-body systems if N is large enough 
(Antoniazzi et al, 2007b; Jain et al, 2007). 

2.3.1 Violent relaxation 

Violent relaxation refers to the initial fast dynamical evolution of collisionless 
systems. It attracted special attention for the following reasons: 

• Contrasting with the Boltzmann equation where the collision term pro- 
vides explicit relaxation, the Vlasov equation does not describe a priori a 
relaxation process. 

• In collisionless systems, it appears that a fast evolution takes place rapidly 
after letting an initial condition evolve, followed by a slower process. This 
phenomenon has been largely discussed for galactic dynamics (Lynden- 
Bell, 1967; Hohl, 1968; Severne and Luwel, 1986; Mineau et al, 1990). 
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We give an example of such a dynamics for the Hamiltonian Mean-Field 
model, considered as a simplification of gravitational dynamics. A molecular 
dynamics in which we follow the magnetization M in the course of time is 
performed (see chapter 4 or 6 for definitions). The initial magnetization is 
Mo = 0.2 and the energy is U = 0.6. In the initial evolution of M, we observe 
strong oscillations, followed by oscillations around 0.3. The run performed with 
N = 1Q 3 particles then drifts towards a higher value. When the number of 
particles is sufficiently high (in our illustration, N = 10 4 ) M remains near the 
value 0.3. It does not experience relaxation apart from the initial evolution 
(0 < t < 50) The fast initial evolution is termed "violent relaxation". In the 
limit where N goes to infinity, this is the only source of relaxation. 



0.45 
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Figure 2.7: An example of violent relaxation in the Hamiltonian Mean-Field 
model. We display the magnetization M in the course of time. After an initial 
value Mq = 0.2, M experiences strong oscillations and stabilizes around a mean 
value of about 0.3 for N = 10 4 while it drifts to a higher value for N = 10 3 . 
The initial regime is termed "violent relaxation" . The simulation is performed 
at energy U = 0.6, ./V is indicated in the legend. 

In order to characterize the state effectively attained by the dynamics after 
the violent relaxation stage and to eventually predict its properties, Lynden- 
Bell proposed a "Statistical mechanics of violent relaxation in stellar systems" 
(Lynden-Bell, 1967). 



2.3.2 Microscopic phase space dynamics 

The Vlasov equation describes an incompressible fluid in phase space. An initial 
condition will be deformed, stretched, and filaments will form but the height 
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and area of the waterbag will remain constant 2 . Already at early times, for the 
Hamiltonian Mean-Field model, we observe a thinning of the contours. Figure 
2.8 displays such a thinning, leading to the formation of filaments. The smallest 
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Figure 2.8: Phase space for the HMF model at time t = 8, with a waterbag 
initial condition of energy U = 0.51 and initial magnetization Mq = 0.05. 

scale eventually attained by the filaments decreases continuously, and it becomes 
impossible to follow that evolution analytically or numerically. 

To clarify the discussion, we here define the meaning of the different scales 
of interest in the study of phase space dynamics. 

1. Microscopic: The microscopic scale refers to the scale describing the fine 
filaments appearing in Vlasov dynamics. It is able to describe fully the 
ideal theoretical evolution of f(6,p). 

2. Mesoscopic: The mesoscopic scale is small with respect to the complete 
phase space, typically too small to be attained by numerical methods but 
larger than the filaments. From the mesoscopic scale, the filaments are 
averaged and not visible. 

3. Macroscopic: The macroscopic scale is the one on which collective phe- 
nomena are observable. For the mean-field models studied in this thesis, 
the characteristic length of interest is typically tt, in relation to the char- 
acteristic scale of the interaction potential. 

The initial condition we consider is a waterbag initial condition (see Fig. 2.4 
for instance), in which / takes only two possible values: or fo- We display 
its evolution through Vlasov dynamics at time t — 8 in Fig. 2.8. In Fig. 2.8, 
filaments are seen to be a prominent characteristic feature in phase space. We 
detail the analysis by displaying the contour lines of the waterbag in Fig. 2.9 

2 The discussion of what happens in numerical simulations is postponed to chapter 4 
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Figure 2.9: Illustration of the different scales of description in phase space. 

A zoom on a particular region is displayed in Fig. 2.10, where we superimpose 
a grid to illustrate the concept of phase space cell. 

The rightmost filament in Fig. 2.10 is too thin to be described by the grid : 
its width is smaller than a grid cell. It is therefore considered microscopic, while 
the leftmost filament is adequately described. 

Figure 2.10 also serves the purpose to illustrate the limitation of the numer- 
ical approach: if we consider the value to assign to grid points of the mesh, we 
can for instance assign to a grid region filled with a filament the area of occupied 
by the filament: for a grid point {9i,Pi), we would get: 

km= / / d9dpf(9,p) (2.20) 

J6i-A9/2 Jp m -Ap/2 

which is a real number between and /q. 

Lyndon-Bell's original idea is the following: the microscopic scale cannot 
be included in the predictive theory, and the prediction will give f(0,p), the 
coarse-grained one-particle probability distribution function, valid at the meso- 
and macroscopic scales. 

2.3.3 Lynden-Bell's entropy 

We describe how to obtain Lynden-Bell's entropy from microscopic phase space 
dynamics. We refer the reader to the following references for more detailed 
explanations (Barre, 2003; Lynden-Bcll, 1967; Campa et al., 2009). 
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Figure 2.10: Zoom of Fig. 2.9. The rightmost filament would need a grid finer 
than the one displayed in order to be correctly described. The scale below our 
"observation" scale is termed microscopic. Lynden-Bell used the analogy with 
the power of a microscope (Lynden-Bell, 1967). 



Lynden-Bell's theory predicts a coarse-grained distribution function /. Whereas 
the microscopic distribution function / can take only two values, or /o, / can 
take any value in that interval. The process is the following: 

1. We decompose phase space in microscopic cells of volume to. 

2. We consider mesoscopic cells containing v microscopic cells. 

3. Considering the "average value" of / in that cell, e.g. 

fj _ n-ifo 



fj" is the microscopic value of / (0 of /o), and rn is the number of micro- 
scopic cells having the value fo- 

4. In that mesoscopic cell, there are ^z^tti ways to place the elements. 

5. There are tj^tt to place the nt's in the mesoscopic cells, where AT = J2i n i 
is the total number of grid cells having the value fo- 
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The result is that the volume W({rii}) occupied by a given configuration is 

TV! -r-r v\ 



w({m}) 



(1/ - ni)V 



(2.21) 



Taking the logarithm, going to the continuum limit and using Stirling's formula 
leads to the Lynden-Bell entropy for a given /: 



/ 



de dp 



f , / 



fo fo 



f 



flnf + 1-f In 1-f 



fo 



f 



fo 



(2.22) 



Maximizing this entropy under constraints (mass, energy, momentum) leads to 
the most probable /. 

Fig. 2.11 illustrates the idea of dividing phase space into cells. f(6,p) in the 
small (microscopic) cells takes only the values or /q. If one considers a group 
of microscopic cells, the average value of f(8,p) is bounded between and /q. 




Figure 2.11: Example of subdivision of phase space into cells. The thick lines 
delimit mesoscopic cells and the small square represent microscopic cells, whose 
value is or fo. From top- left to bottom right, horizontally, the occupation 
levels are: 12, 0, 0, 16, 11, 0, 0, 11 and 13. 
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Chapter 3 



Definition of the 
Hamiltonian Mean-Field 
model 



This chapter is dedicated to the definition of the Hamiltonian Mean-Field (HMF) 
model. We illustrate some of its properties with numerical simulations in the JV-body 
context. 



The Hamiltonian Mean-Field (HMF) model has been introduced in Refs. (In- 
agaki and Konishi, 1993; Ruffo, 1994; Antoni and Ruffo, 1995) [see also (Pichon, 
1994)]. It can be interpreted as either a set of ferromagnetic particles with in- 
finite range coupling or as the simplification of a gravitational model in one 
dimension in which only the first term of the Fourier expansion of the interac- 
tion is kept. The HMF model is meant to give a simplified framework in which 
it is possible to study features of the more realistic systems. 

A phase transition is found for the HMF model in the microcanonical and the 
canonical ensembles, separating homogeneous and non-homogeneous situations 
(Antoni and Ruffo, 1995). Molecular dynamics simulations in the microcanoni- 
cal ensemble show a peculiar behaviour near the phase transition, a feature that 
gave rise to interest in the dynamical behaviour of the HMF model. 

We give in this chapter a definition of the HMF model and illustrate some 
of its dynamical properties. 
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3.1 Definition 



The HMF model is described by the Hamiltonian 

i—l i,j= 1 

where 6i is the position in [— tt; ir[ of the ith particle, pi its momentum, N is the 
number of particles and J is a coupling constant. We consider in this thesis the 
HMF model with J — 1 and do not write this constant anymore. 

The magnetization is commonly used to track the dynamical evolution of 
the system. It is defined as : 



1 N 

M = (M x , My) = — ^(cosfli.sinfli) (3.2) 



and we define M = + . The magnetization also allows one to rewrite 

the Hamiltonian (3.1) in a way explicitly extensive : 



N Pi , N 



^ = Ey + i( 1 - M2 ) • ( 3 - 3 ) 

i=l 

In the following, we make use exclusively of the energy per particle U : 

U=§ . (3.4) 
The equations of motion for the particles are given by : 



W Eli s in (ft 



(3.5) 



= -M(t)s3n(6i-<p(t)), 



where ip is defined so as to satisfy M x + iM y = Me lip (here, i is the complex 
number). Equations 3.5 are similar to the equations of motion for a pendulum 
in which the field intensity and direction are dependent on the time. 



3.2 Dynamical evolution 

We comment on the properties of the TV-body dynamics for the HMF model. 
We focus on the case of waterbag initial conditions that are of relevance in this 
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thesis. In the waterbag initial condition, particles are distributed uniformly in 
phase space according to the following conditions: 



\9i\ < A9 , for all i 
\pi\ < Ap , for all i . 



(3.6) 
(3.7) 



One can characterize the initial distribution equivalently by the energy U = 
^jf- + |(l - M$) and the initial magnetization M = ■ 

Dependence on the initial magnetization M . Statistical mechanics pre- 
dicts for a given energy U a unique value of the magnetization M (Antoni and 
Ruffo, 1995). In contrast, the dynamical evolution of the HMF model depends 
also on the value given to Mq. This is illustrated in Fig. 3.1 where different 
evolutions are observed depending on the value of Mq- A phase diagram sep- 
arating the parameter space (Z7, Mo) in regions with M — and M ^ can 
be computed with the help of Lyndcn-BclPs theory (Antoniazzi ct al., 2007c). 
This is an out- of- equilibrium phase transition, a feature of some systems with 
long-range interaction. 
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Figure 3.1: Evolution of the magnetization M in the HMF model. U — 0.6 and 
TV = 10 4 . Depending on the initial magnetization Mq, we observe significantly 
different values for M. 



Dependence on the number of particles N. The evolution of the mag- 
netization M shows a strong dependence on the number of particles N. After 
strong oscillations in M, the simulations with N = 10 and N = 10 5 particles in 
Fig. 3.2 stabilizes at a value of about 0.35, which is not the value predicted by 
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Figure 3.2: Evolution of the magnetization in the HMF model. U — 0.6 and 
Mo = 0.2 . We observe a qualitatively different evolution depending on the 
number of particles in the simulation. 

statistical mechanics. The simulations with a lower number of particles evolve 
towards a higher value of M. 

The strong dependence of the values attained by M on the initial magneti- 
zation Mq and on the number of particles N represents a conceptual challenge 
for statistical mechanics and nonlinear dynamics. 
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Chapter 4 



Numerical resolution of the 
Vlasov equation 



The nonlinear dynamics resulting from the Vlasov equation is not tractable ana- 
lytically and one needs to resort to a numerical treatment. We introduce the semi- 
Lagrangian method, its numerical properties and its application to the Hamiltonian 
Mean-Field model. We also discuss the case of a model for the Frcc-Elcctron Laser. 



The numerical resolution of the Vlasov equation remains an active research 
domain, despite its foundations being laid as early as in 1976 (Cheng and Knorr, 
1976). This method is a fundamental tool for the understanding of plasma 
physics and is used, for instance, to study the behaviour of tokamaks. The 
major limitations is the resolution (number of grid points) used to describe the 
evolution of the one-particle distribution function (1-PDF). The 1-PDF is an 
infinite dimensional object and its representation on a computer, whatever the 
method chosen, is imperfect. 

The literature on the numerical resolution of the Vlasov equation covers 
mostly the Vlasov-Poisson system of equations describing plasmas in various 
configurations (Shoucri, 2008). Refinements include relativistic descriptions 
( Sonnendruckcr et al., 1999), parallel computations (Crouseilles et al., 2007) 
and the use of wavelets (Gutnic et al., 2004), for instance. 

In this thesis, we use the numerical resolution of the Vlasov equation to 
study one-dimensional models in which peculiar dynamical phenomena have 
been discovered (Dauxois et al., 2002; Campa et al., 2008). The Vlasovian 
approach has played a considerable role in their understanding (Yamaguchi 
et al., 2004; Antoniazzi et al., 2007c; Antoniazzi et al., 2007a); meanwhile, the 
performed simulations have remained almost only of the molecular dynamics 
type. The study of direct Vlasov simulations is an essential addition to this 
approach and requests a detailed analysis of their properties. 
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Conversely, taking profit of the relative simplicity of our mean-field mod- 
els and of the growing interest for their Vlasov dynamics in the literature, we 
compare our simulation to recent studies on the Vlasov- Poisson system that are 
concerned with the limitation of the numerical scheme. These comparisons aim 
at making the Hamiltonian Mean-Field model a test-case for Vlasov algorithms. 

This chapter is structured as follows: we review the problem of the numerical 
resolution of the Vlasov equation in section 4.1. We then proceed to detail the 
resolution of the ID advection equation that forms the basis of the complete 
algorithm in section 4.2 We present the algorithm for the complete ID Vlasov 
problem, following the work of Cheng and Knorr (Cheng and Knorr, 1976) 
and Sonnendrucker et al ( Sonncndrucker et al., 1999) in section 4.3 and detail 
its application to the Hamiltonian Mean-Field model (section 4.4) and to the 
Free- Electron Laser (section 4.5 ). 

4.1 Overview of the problem 

In order to study the evolution of f(0,p;t) , we need to choose a method to 
represent and store this distribution function on a computer. Several points of 
view exist, which we briefly recall: 

1. Lagrangian Methods : 

/ is here computed with the knowledge of a set of particles. These particles 
represent themselves a number of microscopic elements and are evolved at 
each time step. This methodology is called particle- in-cell (PIC). 

2. Eulerian Methods : 

The computation is based on / directly. The function is stored on the 
computer, and serves as the basis for the computation of the force and 
of the macroscopic quantities. A number of methods exist to store / 
(regular Cartesian grid, irregular grid, a Fourier basis, wavelets), and the 
implementation of the algorithm is tailored to the storage method. 

The main advantage of the PIC method is its ability to treat the full 6D 
problem in terms of storage capacity [see for instance Ref. (Hockney and East- 
wood, 1988)]. Its major drawback is to suffer from an important noise level 
arising from the sampling in phase space. Another fundamental difference is 
that, in the case of Eulerian methods, one needs to run a single simulation by 
physical situation. Since there is no random generation of particles correspond- 
ing to the macroscopic situation, the averaging over several realisations is not 
needed. 

Given the modest computational cost in storage as well as in CPU-time 
for ID systems, we can afford the choice of an Eulerian method. The highest 
resolution that we will use, 8192 x 8192, requires w 540MB of memory for a 
copy of / in double precision, making the computation accessible, memory-wise, 
to single-CPU computers. The use, in our algorithm of /, a copy of / and the 

lr The notation and description of physical quantities is discussed in chapter 2 
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storage of the second derivative needed by the spline interpolation elevates the 
memory cost by a factor of three. 

In order to assess the quality of the numerical resolution, we make use of Li 
norms, which are ideally conserved. These conserved quantities are not expected 
to hold exactly (up to machine precision) in the numerical resolution. We will 
compute the following quantities: 

• Li[f] corresponds to the normalization of /. 

• Deviation of Li [/] indicates the presence of numerical dissipation. 

• The energy U[f] is also a constant. 

The limitation of the resolution to describe / during the development of 
filaments comes from the fact that the ideal / needs a microscopic phase space 
definition (see section 2.3.2 for a definition), which is impossible. The problem 
of filamcntation is detailed in section 2.2.5 

A rough criterion based on the free-streaming situation gives an estimation of 
the time during which the numerical simulation is able to follow the development 
of the filaments. We recall the steps given by Canosa et al (Canosa et al., 1974). 
The explicit solution to the free-streaming Vlasov equation is written in terms 
of a Fourier decomposition, for an initial condition perturbed at wavenumber k: 

f(6,p;t) =eiq>(ik6)exp(-ikpt) . (4.1) 

Equation (4.1), specialized for a discrete set of momenta p m = mAp, reads : 

f(9,mAp;t)=cxp(ik9)exp(-ikmApt) . (4.2) 

Equation (4.1) indicates that for any given to, the distribution function comes 
back to its initial value at a time t m — k ^ Ap - Considering the resolution for a 
set of values to, one finds that after a time : 

r " = i; • <«> 

the complete distribution function / experiences a recurrence, clearly indicating 
a limitation of the resolution procedure. Tr is called the recurrence time. 



4.2 The advection equation 

To solve the Vlasov equation, we make use of the time-split algorithm intro- 
duced in the work of Cheng and Knorr (Cheng and Knorr, 1976). Only a 
one-dimensional advection problem needs to be solved in this situation. 
The ID advection with a constant velocity field reads: 

at ox 
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Considering the problem on a fixed grid, a simple way to solve numerically this 
problem is a finite-difference explicit scheme. We note ff the value of f(xi \ sAt) 
where Xi = x + iAx, i e [1; N x ], and At is the time step. 
The finite-difference version of Eq. (4.4) is: 

fs+l _ fs f s _ p 

Jl - Jl +c Jl = 0, (4.5) 

At Ax 

giving an explicit formula for // + . 

This formulation allows for a simple computation of / at time s + 1 given / 
at time s. It is however strongly dissipative. Dissipation in this situation refers 
to an artifact of the algorithm acting as a diffusion operator. A simple finite- 
difference scheme as this one can be used only in situation where a physical 
diffusion operator is present, as in Fourier's heat equation. Another important 
property of this scheme is the bound on the value of the time step, the so-called 
Courant-Friedrichs-Lewy condition [see for instance Ref. (Zwillinger, 1998)] 

At , , 

c-<l . (4.6) 

Let us introduce the method of the characteristics that will be used further 
on. This algorithm is based on the fact that, for Eq. (4.4), the value of f(x; t + 
At) correspond to the value of f(x*;t) where x* is located at the foot of the 
characteristic curve. We focus on the ID problem, but the presentation is valid 
in higher dimensions with the same notations. 

x* is computed (Zwillinger, 1998) according to : 

/t-At 
dt c = x-cAt (4.7) 

allowing to write an explicit solution to Eq. (4.4) on a grid as 

ft +1 =f s {x t -cAt) (4.8) 

In general, x* does not lie on the storage mesh, meaning that we need an 
interpolated value of / to proceed. The choice of the interpolation method has a 
great impact on the numerical properties of the resolution, as detailed in several 
articles (Filbet et al., 2001; Arber and Vann, 2002). 

In the rest of this work, we use the cubic spline method (Press et al., 1996) 
that gives a good compromise between numerical dissipation, computational 
time and complexity. 



4.3 The semi-Lagrangian method for the Vlasov 
equation 

Solving the Vlasov equation in ID amounts to solve a 2D partial differential 
equation, one for the position and one for the velocity. Cheng and Knorr have 
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f*{6i, Pm ) = f s (6i - Pm At/2, Pm ) 



already proposed in 1976 a simplification via a time-splitting algorithm (Cheng 
and Knorr, 1976). Their work laid the foundation for further works. 

Cheng and Knorr introduced the following discretization, accurate to second 
order: 

f +1 (8,p) = f s (0-At(p-^F*(6)At),p-AtF*(9)) (4.9) 

where 9 = d — pAt/2 and F* is the force field computed at half a time step. We 
use the notation f s {9 1 p) — f(9,p; s At). 

A practical implementation to compute Eq. (4.9) on a numerical mesh is 
the following : 

1. Advection in the 
^-direction, 1/2 time step 

2. Computation of the force field 
for f* 

3. Advection in the 
p-direction, 1 time step 

4. Advection in the 
^-direction, 1/2 time step 

Each of these steps is performed for the whole grid before going to the next 
one. A sketch of these steps is given in Fig. 4.1. This method, involving a fixed 
mesh and trajectories along the characteristics backwards in time is said to be 
semi-Lagrangian. It is used in fluid dynamics, e.g. for weather forecasts (see 
(Staniforth and Cote, 1991)). 



r*(6i,p m ) = f*{6i,p m - F*{6i)At) 
f s+1 {0 t , Pm ) = f**(6i-p m At/2,p m ) 



a) 




Figure 4.1: The semi-Lagrangian algorithm: a) illustrates the computation of 
the value of / at the grid points indicated by the dots from the feet of the 
characteristic curves (the start of the arrow) . b) and c) represent steps 1 and 3 
of the time-split algorithm. 

In the framework of the Vlasov equation, especially the Poisson-Vlasov sys- 
tem, it was introduced by Sonnendriicker et al ( Sonnendrucker et al., 1999), 
also associated with cubic splines. Let us mention that the semi-Lagrangian 
method is symplectic if one neglects the approximation present in the inter- 
polation step. Watanabe and Sugama proposed a generalization of symplectic 
algorithms for the Vlasov equation and presented the result up to sixth order 
in time in Ref. (Watanabe and Sugama, 2005). 
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4.4 The study of the Hamiltonian Mean-Field 
model 



We discuss in this section the computations related to the properties of the 
simulation for an initial Gaussian homogeneous distribution, then we check the 
stability criterion of Ref. (Yamaguchi et al., 2004). References (Galeotti and 
Califano, 2005; Califano et al., 2006) propose original analyses for the numerical 
resolution for the Vlasov-Poisson system of equations of relevance for plasma 
physics. We comment our results in the same spirit. 

4.4.1 The Vlasov equation for the HMF model 

The Vlasov equation for the Hamiltonian Mean-Field (HMF) model is : 

?L + M dv if]df _ 

dt P d8 d9 dp 
V[f}{6) = 1 - M x [/] cos 6 — My [/] sin 6 
M x [f] = / dOdp f cos 6, 

M x [f] = jdBdpfwa.8 . (4.10) 
Equations (4.10) preserve the total energy U : 

U(t)[f} = Jdedpf(0,p;t)(^ + ^(l-M x [f}co S e-M y [f}smO)j (4.11) 
and the total momentum P : 

P(t)[f] = f d6dpf(6,p;t)(p). (4.12) 



4.4.2 Comparison of Vlasov and molecular dynamics 

We compare the dynamics of the Vlasov resolution with TV-body simulations. 
The Vlasov equation is a good description of the system up to a validity time 
that depends on the number of particles. Braun & Hepp (Braun and Hepp, 
1977) prove the convergence to the Vlasov equation for weakly coupled systems. 
Jain et al (Jain et al., 2007) detail the validity time in the case of the HMF for 
homogeneous initial conditions and quasi-stationary states. 

In this paragraph, we compare a Vlasov simulation with a sligthly perturbed 
initial condition (4.13) with V-body simulations setup in the same fashion, but 
with no perturbation. The instability is triggered by fluctuations of order 1 / y/~N 
in the magnetization. The energy is set to U = 0.51. For the Vlasov simulation, 
the recurrence time is Tr w 357, which is longer than the simulation. 

The initial dynamics is similar for all simulations, but after t w 30, the run 
with N = 1000 particles clearly separates from the two other runs, confirming 
that the Vlasov dynamics corresponds to a high number of particles. 
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Figure 4.2: Comparison of the magne- 
tization M in Vlasov and molecular dy- 
namics simulations for the HMF model. 
Run G512 is described in section 4.4.3, 
the energy is U — 0.51. 
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t 



Figure 4.3: Zoom on Fig. 4.2, same leg- 
end. The simulation with N = 10 6 
stays close to the Vlasov simulation 
for a longer time than the one with 
N = 10 3 . 



Let us remark that the use of the energy per particle U — H/N and of the 
average magnetization M = £\ (cos0j,sin#j) allow to cast the results for 
various values of TV in a simple manner. 

4.4.3 Properties of the numerical solution 

We run simulations starting with a Gaussian profile, slightly perturbed : 



/(M = 




(4.13) 



at an energy of £7 = 0.51, e= 10 4 . j5 is computed from U, for an homogeneous 
state : /? = \ u\/2 ■ At U — 0.51, the Gaussian profile is unstable; M grows, 
then oscillates around a mean value to which it eventually relaxes. 

The following parameters are used : At = 0.1, size of the box in the p- 
direction [—4.5 : 4.5]. The grid size is varied from N$ — N p = 64 (G64) to 
N e = N p = 512 (G512); we call these runs G64, G128, G256 and G512. 

We detail in Table 4.1 the conservation properties of the algorithm. The 
conservation of U could be improved by decreasing At, implying more compu- 
tational time. This is the sole constraint on At as the semi-Lagrangian method 
has no Courant condition. 

M displays a similar behaviour for all values of the resolution until t ~ 60, 
except for the lowest resolution run G64; M has been checked to correspond to 
JV-body dynamics with N = 10 6 for short times (data shown in section 4.4.2). 
M grows up to a saturation value, identical for all runs, then oscillations take 
place and eventually damp to a stable value. 

The damping of the oscillations is strongest in G64, they disappear at t s» 70. 
This damping is present in the other runs as well, a phenomenon that diminishes 
when the grid size is increased. This fact is to be put in correlation with the 
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run 


Li(t)-Li(0) 
maX L l( 0) 


mnv U(t)-U(0) 

max U{Q) 


G64 


1.6 10~ 4 


3.2 1(T 4 


G128 


2.1 10~ 5 


2.0 10~ 4 


G256 


7.2 10~ 6 


2.0 10~ 4 


G512 


2.8 1CT 6 


2.0 10~ 4 



Table 4.1: Conservation properties for different grid sizes. Increasing the grid 
size allows for a better conservation of L\ . The conservation of U depends also 
on the time step (second order scheme) and could be improved by decreasing 
At. 




400 



Figure 4.4: Evolution of M for the Gaussian initial condition given by Eq. (4.13) 
(U = 0.51,6 = 10~ 4 ). A magnification of the same curve is presented in the inset. 



evolution of L2 (Fig. 4.5) : The decrease of L2 means that the small scales are 
not well described anymore by the numerical procedure. This decrease is the 
strongest in G64 where L2 does not reach a plateau and diminishes continuously. 
The shape of /, a vortex rotating in phase space, is giving the oscillations of 
M. When this structure is smoothed out by the lost of the small scales, see 
Fig. 4.10, M tends to a constant value. The run G512 displays a filamentary 
structure and oscillations in M up to the end of the simulation. 

The absence of oscillations in M are directly linked to the fact that / becomes 
invariant under rotations in phase space in the process of numerical dissipation. 
Macroscopic quantities (such as M) cannot vary in time anymore, except for 
the small effects of the ongoing dissipation (see the continuous decrease of L2 
for run G64 in Fig. 4.5). This phenomenon may differ significantly from the 
ideal evolution of / and M : in a number of situations, oscillations persist due 
to the presence of coherent structures; this is only observed in run G512. 

These observations suggest a limit in time of the validity of the simulations. 
In the case of the free streaming, the evolution of / leads to a recurrence after 
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Figure 4.5: Evolution of L2, same runs as Fig. 4.4 



a time Tr — 27r/(fcAv), called adequately "recurrence time" (k is the mode 
number of the initial perturbation and Av is the grid spacing in velocity space, 
see section) at which / recovers its initial value (Canosa et al., 1974). In situa- 
tions close to free streaming, for instance in the linear Landau damping, Tr is 
clearly identified in the evolution of the first mode of the electric field (Canosa 
et al., 1974). In other situations, one can only use Tr as a rough estimate and 
cannot identify its effects as sharply. It however agrees with observations made 
otherwise observing M(t) or cuts in /. For the run G64 the recurrence time Tr 
is approximately 45, discarding the simulation in Fig. 4.8 where t — 64, agreeing 
with the observed discrepancy; for the run G512, we are close to the simulation 
time : Tr 357. 

Asymptotically, all runs tend to a similar value of M , which we characterize 
with the standard deviation of M, ctm, and the mean value across simulations, 
M : <j m /M = 1.1%, or 0.7% if we leave G64 aside. 



4.4.4 Stability of the initial condition 

The stability of homogeneous initial condition has been studied in Ref. (Ya- 
maguchi et al., 2004). The authors of Ref. (Yamaguchi et al., 2004) devise a 
criterion of formal stability for any homogeneous initial condition. We consider 
a distribution function /o of the form : 

/o(fl,p) = ^o(p) ■ (4.14) 

Given a fixed point of the form given in Eq. (4.14) of a conserved functional 
F[f], 

F[f] = C s [f] - PHvif] -»[ f(0,p) d0 dp , (4.15) 

where C s is a Casimir (depending on a concave function s) and Hy the energy 
functional, one computes its variation AF[5f] — F[fo + Sf) — F[fo] . The sign of 
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the second order term of AF determines the formal stability of fa- s disappears 
from the computation and criterion, depending only of /o, becomes : 



I[fo 



1 + 7T 



+ °° d p f o (0,p) 

oo P 



dp > /o is formally stable. 



(4.16) 



The criterion defines, for a given family of /oi a critical energy U* below which 
the initial condition is unstable. 

We perform a check for the waterbag and Gaussian initial profiles. The 
waterbag is the following profile, for a given U: 



0, else. 



(4.17) 



where we have again applied a small perturbation of order e to trigger the 
instability, e is set to 10 -4 for all runs. 

Starting from both side of the critical energy (U* = 7/12 for the waterbag 
and U* = 3/4 for the Gaussian), we observe in Figs. 4.6 and 4.7 the change in 
stability. For values of U above U*, M stays at the same order of magnitude or 
decreases. 




10 20 30 40 50 60 70 
t 
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Figure 4.6: M(t) for initial Gaussian profiles. M grows if U is below U* = 3/4. 



In the case of the waterbag, Yamaguchi et al (Yamaguchi et al., 2004) 
compute the exponential growth rate of M(t). We compare that prediction 
with the simulation results of Fig. 4.7 [fitting logM(i)] and found a very 
good agreement. The exponential growth rate depends on the energy U : 
A = y^6 (U* — U). The theoretical and numerically computed values are given 
in couples (A+u,Anum) f° r U = 0.52,0.55 and 0.57 respectively : (0.61,0.61), 
(0.43,0.43) and (0.32,0.31). 
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t 

Figure 4.7: M(t) for initial waterbag profiles. M grows if U is below U* — 
7/12 ss 0.58. The slope before saturation of the three first runs is used to 
compute the exponential growth factor A. 

4.4.5 Effect of the numerical grid size 

The results of section 4.4.3 have considered different grid sizes for the same 
problem; f(9,p) was inspected visually, but most of the analysis focused on 
macroscopic quantities : L\, L2, U and M. Galeotti, Califano and Mangeney 
(Galeotti and Califano, 2005; Califano et al., 2006) compared cuts in /, i.e. 
f(6,p*) where p* is fixed or f(9*,p) where 9* is fixed, after smoothing of all 
runs to the lowest accuracy. This method explores with more precision the 
structure of /. For runs G64 to G512, we display f(9,0) in Figs. 4.8 and 4.9 
at two times : t — 64 when runs G128 to G512 still have a similar M, and at 
t = 400 in order to discuss the asymptotic evolution. 

At t = 64, although the values for Li are different, which may imply a 
different evolution of /, we find peaks in f(9,0) at identical locations. This 
indicates that / behaves similarly. Indeed, the force on the particles (or phase 
space elements) only depends on M. Thus, as long as M is identical for the 
runs, so is the force term. What this means is that an equal M(t) gives a 
similar f(9,p;t), up to the accuracy of a given grid on which / is represented; 
this claim is supported by Fig. 4.8 where the runs G128-G512 agree very well, 
in correlation with the agreement of M(t) shown in Fig. 4.4 up to i « 64. 

At t = 400, G64 is completely smooth, and the other runs still show a small 
filamentary structure. The location of the maximum is the same, a fact that 
contrasts with most results of Ref. ( Califano et al. , 2006 ) for a Vlasov plasma, 
and that indicates a similar asymptotic state. Regarding the position of the 
final vortex, the initial dynamics of the runs G64-G512 forms only one vortex, 
whose location is stable, immediately after the growth of the initial perturbation 
of the initial condition (4.13); an identical position is found in Ref. (Califano 
et al., 2006) in a particular case (final runs D,E and F). On the other hand, the 
agreement in the cut of / for the runs G128-G512 in Fig. 4.8 is better in our 
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Figure 4.8: Cut in the function / for p* — at time t = 64. Same runs as in 
Fig. 4.4. The locations of the peaks is similar in all runs. We observe that the 
heights of the peaks in run G64 stand apart. 
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Figure 4.9: Cut in the function / for p* = at time t — 400. Same runs as 
in Fig. 4.4. The location of the maximum is the same in this asymptotic state, 
and the value of the maximum is very similar for all runs. 

case and reflects the similarity in M(£), as discussed in the previous paragraph. 

We learn from these observations, in the case of the HMF model, the fact 
that one observable, M , describes not only the macroscopic evolution, but also 
monitors the "quality" of / (i.e. the macroscopic evolution and the fact that 
/ has conserved a shape which is not invariant under rotation in phase space 2 ; 
We point out that M is not meant to replace the more complete information 
contained in /). The quantity M is similar to the energy spectrum for a Vlasov 
plasma, the energy spectrum of the HMF model containing only one mode. The 

2 See discussion on the oscillations in paragraph 4.4.3 
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development of Vlasov codes may benefit from the simplicity of using the single 
quantity M . This is especially useful in the process of comparison between 
different algorithm that is often carried out in numerical analysis (for Vlasov 
plasmas, see for instance Ref. (Arber and Vann, 2002)). 




-3-2-10123 -3-2-10123 




3-2-10 1 2 3 




Figure 4.10: Phase space contour lines : a) and b) depict the run G64 at times 
32 and 64. c) and d) depict the from run G512 at the same times. While a) 
still presents the correct shape as compared to c) , b) has lost most of the spiral 
shape present in d). 



4.4.6 Conclusions 

We presented in this chapter a detailed study of Vlasov numerical simulations 
for the Hamiltonian Mean-Field (HMF) model, that have been compared to 
previous work on the HMF model and to Vlasov results in plasma physics. We 
gave evidence that the mean-field order parameter M of the HMF model is 
sufficient to compare different runs in terms of phase space evolution as well as 
dissipative (i.e. damping) properties. 

Illustration was given of the "non- Vlasov" limit of the numerical simulation 
which we need to keep in mind to analyze simulations results carefully, as stated 
recently by other authors (Galeotti and Califano, 2005; Califano et al., 2006; 
Antoniazzi et al., 2007a). 

Finally, we hope that cross research between long-range systems obeying a 
Vlasov equation and Vlasov plasmas may benefit to both fields by providing the 
first with a new tool and the latter with new insight on the numerical method. 
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4.5 Application to the Free-Electron Laser 



The Vlasov equation for the Colson-Bonifacio model for the Free-Electron Laser 
(FEL) reads (see chapter 5): 

df df df 

f_ = -p-± + 2(A x cos6-A v smO)-±, (4.18) 

8 A r 

-^r = d0dpf cos 9, (4.19) 

8 A f 

= -JdBdpfwa.9. (4.20) 

These equations are very similar to the ones of the HMF model, but the mean- 
field term is not computed self-consistently. A x and A y are evolved through 
Eqs. (4.19) and (4.20). An initial condition for / is accompanied by an initial 
value for A x and A y . When the initial / is homogeneous, almost vanishing value 
for A x (0) and A y (0) correspond to self-amplified spontaneous emission (SASE). 
We will set A x (0) = 10~ 4 and A y (0) = in the following simulations, in order 
to trigger the instability. 

The quantities b x [/] and b y [/] allow to compute the bunching of the electrons 
b=Jbl + 62; 

b x [f}= J d6dp /cos0 ; b y [f] = J dOdp fsmO (4.21) 



We do not perform an in-depth study as for the HMF model, but present 
the general properties of simulations of interest for chapter 5 . 



4.5.1 Time splitting 

In order to integrate the system (4.18), (4.19), (4.20), a time splitting needs to 
be defined. After testing several possibilities, we have settled for the following: 

1. Advection in the ^-direction, 1/2 time step 



2. Computation of b* x and b* 

3. Advection in the p-direction, 1 time step 

4. Advection in the ^-direction, 1/2 time step 

5. Computation of A s + l = A% + \(V X + b x ) and A s y +1 =A y -\{b* y + b y ) 

This time splitting scheme avoids energy drifts, which is an important property 
for the situations encountered in paragraph 5.6 
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4.5.2 Properties 

Figure 4.11 displays the evolution of / = A\ + Ay as a function of z for a 
waterbag initial condition with 6 = 0.05 and energy U = 0.10, for grid sizes 
going from N e = N p = 64 (WB64) to N e = N p = 512 (WB512). As in the case 
of the HMF model, the initial dynamics is similar, and we observe a growing 
disagreement between the runs. The average final intensity is similar for all grid 
sizes. 




Figure 4.11: Evolution of I = A\ + Ay as a function of z in Vlasov simulations 
for the Free-Electron Laser. 

The oscillations in I(z) follow a similar ordering as for the HMF model, 
indicating that the run with the largest grid size run is able to describe details 
that are too fine for the other runs. The added precision makes a qualitative 
difference in the behaviour of I. 

Table 4.2 presents the conservation properties of the simulations and Fig. 4.12 
displays the evolution of the Li norm. The normalization is well preserved and 
the energy conservation does not improve significantly with the resolution. 
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run 


Li(t)-Li(0) 
maX L l( 0) 


U(t)-U(0) 

max uw 


WB64 


1.2 10~ 3 


2.3 10~ 2 


WB128 


6.8 1CT 5 


3.8 10~ 4 


WB256 


1.0 10~ 5 


3.9 10~ 4 


WB512 


4.9 10~ 6 


4.3 10~ 4 



Table 4.2: Conservation properties for different grid sizes. Increasing the grid 
size allows for a better conservation of L\ . The conservation of U depends also 
on the time step At and could be improved by decreasing At. 
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Chapter 5 



Out-of-equilibrium 
dynamics in the 
Colson-Bonifacio model for 
the free-electron laser 



The Free-Electron Laser is a device in which a coherent beam of light is produced 
from a beam of electrons. A simple model for the wave-particle interaction describes 
the Free-Electron Laser and manifests a complex nonlinear dynamics. The Vlasovian 
formulation is given and we perform the corresponding simulations. In the context 
of the Free-Electron Laser, we describe the out-of-equilibrium phase transition taking 
place between lasing and non-lasing phases. We interpret the transition with the help 
of Lynden-Bell's theory and comment on the structures appearing in phase space. 



Mean-field models have been widely studied as paradigmatic representatives 
of the important class of systems subject to long-range coupling. In the simplest 
scenario, N particles are made to interact in one dimension, subject to a varying 
field which is self-consistently sensitive to the individual trajectories. A global 
network of connections is hence driving the dynamics of every constituting el- 
ement, as it certainly happens for more realistic settings where e.g. gravity or 
unscreened Coulomb interactions are at play. 

Free-electron lasers (FELs) are lasing devices consisting of a relativistic beam 
of charged particles, interacting with a co-propagating electromagnetic wave. 
The interaction is assisted by the static and periodic magnetic field generated 
by an undulator. FELs admit a mean-field description in term of the so-called 
Colson-Bonifacio model (Colson, 1976; Bonifacio et al., 1987), which captures 
the essence of the collective wave-particle dynamics. 
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The FEL displays quasi-stationary states (QSS) (Barre et al., 2004) and 
is known to possess a threshold in energy above which the intensity drops to 
zero. These QSSs bear an extraordinary conceptual importance as they poten- 
tially corresponds to the solely experimentally accessible regimes, for instance 
in plasma physics. 

The QSS in the Hamiltonian Mean-Field model (HMF) have been explained 
by resorting to a maximum entropy principle, pioneered by Lynden-Bell in as- 
trophysical context, and fully justified from first principles (Lynden-Bell, 1967). 

The Lynden-Bell protocol, also termed violent relaxation theory, has already 
proved effective in predicting the saturate intensity of a FEL (Curbis et al., 
2007; Barre et al., 2004). However, no detailed study has been carried out for 
the FEL case aiming at unravelling the possible existence of out-of-equilibrium 
transitions of the type mentioned above. Are these transitions ubiquitous in 
mean field dynamics and, in this case, can we provide a consistent interpre- 
tative framework for their emergence? This chapter is dedicated to answering 
such questions, and makes reference to the specific FEL setting. We also stress 
that the Colson-Bonifacio model of FEL dynamics can be regarded as a gen- 
eral formulation for all those applications where the complex interplay between 
particles and fields is well known to be central, e.g. electrostatic instabilities in 
plasma physics (Elskens and Escande, 2003). As a closing remark, we notice 
that, on the practical implication side, by disposing of reliable predictive tools 
on the system evolution, one can aim at guiding the system towards different 
experimental regimes. 

This chapter is structured as follows: Section 5.1 casts the model and its 
Vlasov formulation, whose dynamics is studied in section 5.2. The Lynden-Bell 
theory is formulated in section 5.3 and forms the basis of the analysis performed 
in section 5.4 and summarized in section 5.5. The structures appearing in phase 
space are discussed in more detail in section 5.6 



5.1 The FEL model 



The Colson-Bonifacio model for the FEL dynamics describes the coupled evo- 
lution of the electrons with a co-propagating wave (Bonifacio et al., 1990). The 
equations read: 




< ^ = -Ae^-A'e-* 9 *, (5.1) 

\ dz N e ' 

where 9j stands for the particle phase with respect to that of the optical wave, Pj 
being its conjugate normalized momentum. The complex quantity A = A x +iA y 
represents the transverse field and N the number of electrons composing the 
electron bunch. We assume here that electrons are perfectly resonant with the 
ponderomotive field generated by the optical wave and by the undulator. Should 



64 



that condition not be met, an additional term quantifying the deviation from 
this ideal behaviour can be added to the above description. We make use of 
that possibility in sections 5.6.2 and 5.6.3. 

In Eqs. (5.1), z labels the longitudinal position along the undulator and it 
effectively plays the role of time. The intensity of the laser field is / = A 2 X + Ay. 
As it can be seen from the last of Eqs. (5.1), the bunching term, b = e~ lBi , 

is the source of wave amplification. The bunching quantifies the degree of lo- 
calization of the electrons in the generalized space of their associated phases. 
The above discrete system of equations admits a Hamiltonian formulation to 
which we shall make reference as the TV-body model. In the TV — > oo limit, 
the system (5.1) converges to the following Vlasov-wave set of equations (Barre 
et al., 2004): 



dA x 

dz 
8Ay 

dz 



d9dp / cos 6*, 

J d9dp / sin 9. (5.2) 



Eqs. (5.2 ) can be simulated numerically, allowing us to monitor the evolution 
of the phase space distribution function f(0,p) along the z axis. The numerical 
algorithm is given in chapter 4, with details regarding the FEL in section 4.5. 

The results of the numerical integration are also checked versus TV-body sim- 
ulations and shown to return a perfect matching on relatively short time scale, 
for large enough values of TV. On longer times, finite-TV corrections do matter. 
The discrete system is in turn sensitive to intrinsic granularity effects, stemming 
from the intimate finiteness of the simulated medium, and progressively migrate 
from the Vlasov state towards the deputed equilibrium configuration. When in- 
creasing its size, the system spends progressively more time in the Vlasov-likc, 
out-of-equilibrium regime. Formally, in the TV — > oo limit, it never reaches 
equilibrium, being permanently trapped in the QSS. 

As previously anticipated, our study is hence ultimately concerned with 
the emergence of QSSs, in a context where particles and waves evolve self- 
consistently. We shall be particularly interested in elucidating the occurrence of 
out-of-equilibrium phase transitions via dedicated numerical simulations, and in 
substantiating our claims analytically. This study allows us to virtually extend 
the conclusion of Ref. ( Antoniazzi et al., 2007c) to a broad spectrum of poten- 
tially relevant applications, beyond the specific case under inspection. Among 
other, it is again worth mentioning plasma physics: A formulation equivalent 
to model (5.2) is in fact often invoked, when studying the collective effects of 
beam-plasma dynamics (Elskens and Escande, 2003). 
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5.2 On the initial conditions and their subse- 
quent dynamical evolution 



Let us turn to discussing our results, as obtained via numerical integration of 
(5.2). In order to make contact with the investigations reported in (Antoniazzi 
et al., 2007c), we shall employ in the following a two-dimensional water-bag 
initial condition in phase space, which can be seen as a rough approximation of 
a smooth Gaussian profile. A (rectangular) waterbag is formally parametrized 
by two quantities, namely the semi- width of the spanned interval in phase, AO, 
and its analogue in the momentum direction, Ap. The corresponding expression 
for / can be cast in the form (see also Fig. 5.5 top- left): 

r /o if bi < a p , 

/(M= \0\<A6, (5.3) 

[ otherwise. 

The initial conditions can be also characterized by defining 

{l _ sin A8 
b °Z ^T' (5-4) 
£ 6 ' 

where bo is the initial bunching, and e the initial average kinetic energy per 
particle. Notice that we access all possible values of the bunching b e [0, 1] by 
properly tuning AO, and all positive energies e by varying Ap. Here, we limit 
our discussion to the case of vanishing initial optical field, Iq ~ 0, the relevant 
parameter space being therefore solely bound to the plane (&o, e )- 

This choice of initial parameters allows for a rich out-of-equilibrium dy- 
namics. We insist nonetheless on the fact that we are specializing on a given 
bi-dimensional subset, and deliberately ignore the third, in principle available, 
direction of the reference parameter space. Quantifying the role of such an addi- 
tional degree of freedom ultimately amounts to investigate the so-called seeded 
configuration (Doyuran et al., 2001; Dc Ninno et al., 2008) and will be the 
subject of future work. 

Let us start by discussing the simplest scenario, where the initial beam of 
particles is uniformly distributed over [— ir;n]. From a physical point of view, 
this amounts to specialize to the case of Self- Amplified Spontaneous Emission 
(SASE, (Bonifacio et al., 1984; Brinkmann, 2006)), where &o = and no seed 
is applied externally. Such a choice was also considered by Barre et al. (Barre 
et al., 2004) and Curbis et al. (Curbis et al., 2007), where the dependence of 
the system evolution on the energy was numerically monitored within the N- 
body discrete viewpoint. Interestingly, 6o = is a stationary solution of the 
Vlasov system: a local perturbative calculation can hence be straightforwardly 
implemented so as to investigate its inherent stability (Bonifacio et al., 1990); 
the calculations are detailed in appendix A. For e < 0.315, an instability occurs: 
both the wave intensity and the bunching factor rapidly grow, before relaxing 
towards an oscillating plateau. The average value of / reached in the oscillating 
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regime is called the saturated intensity I. This behaviour is displayed in Fig. 
5.1, where the simulations with e > 0.315 do not show an amplification of /. 
This is a well-known property, indeed correctly reproduced by our numerical 
simulations, and which first signals the existence of phase transitions, of the 
type depicted in (Antoniazzi et al., 2007c). To further corroborate our guess on 
the bo = behaviour, we turn to measuring the saturated intensity / as function 
of the energy e, where / stands for the mean of / after saturation. It is here 
computed during four oscillations of I. 

As shown in Fig. 5.2, I rapidly shrinks, when increasing the energy e, until a 
critical value is reached where a sudden transition to I ~ is observed, bearing 
the characteristic of a first order phase transition. This is a further point of 
contact with the analysis carried out in (Antoniazzi et al., 2007c) for the HMF 
toy model. 




Figure 5.1: (Color online) As a result of the numerical integration of Eqs. (5.2), 
the evolution of / as a function of z is reported for different choices of the energy 
e (e = 0.05,0.1,0.2,0.3,0.4); 6o = 0. Symbols pinpoint the position of the first 
peak in the intensity time series, thus returning an indication on the saturation 
time. Notice that for e > e c = 0.315, the peak is found for / 1 : the 
corresponding initial conditions are hence stable, and no instability develops. 

Motivated by these findings, and to push the analogy with the HMF setting, 
we consider bunched initial distributions. From a physical point of view, this 
choice is relevant to the case of FELs working in the so-called harmonic genera- 
tion regime (Doyuran et al., 2001; De Ninno et al., 2008). The positions of the 
particles are here initially assigned in order to uniformly span a limited portion 
of the allowed support, symmetric with respect to the origin, controlling the 
associated bunching via Eq. (5.4). 

The inhomogeneous (bo ^ 0) distribution in phase space is by nature non- 
stationary. The non-zero bunching amplifies the wave that rapidly acquires a 
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Figure 5.2: (Color online) Saturated intensity, /, vs. e for different choices of 
the initial bunching b = 0.0, 0.05, 0.20, 0.50 and 0.90. 




e 

Figure 5.3: (Color online) Same as Fig. 5.2, for the saturated bunching, b. 

non-zero value rapidly and, thanks to its initial velocity profile, the waterbag 
changes its shape with no possibility to remain stationary. The Vlasov dynam- 
ics can however smooth it out into a homogeneous distribution (b = 0, I = 0), 
possibly not of the waterbag type, or evolve to a bunched situation. The sat- 
urated mean-field average intensity I vs. the energy parameter e is depicted 
in Fig. 5.2, showing the newly collected data for different values of &o > to 
the reference profile relative to bg = 0. In all cases the intensity is shown to 
decrease, as the energy increases. Importantly, for small values of bo, an abrupt 
transition is observed, which can be naively interpreted as of the first-order 
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type. For larger values of bo, the observed transition becomes smoother, such as 
for a second-order one. A substantially identical scenario holds for the bunch- 
ing, which evolves towards an asymptotic plateau b, also sensitive to the e and 
60 parameters, see Fig. 5.3. This scenario points towards a unifying picture 
on the emergence of out-of-equilibrium phase transitions within the considered 
class of mean-field Hamiltonian model. As previously anticipated, the Lynden- 
Bell theory of violent relaxation was successfully applied to the HMF problem, 
allowing one to gain a comprehensive understanding on the out-of-equilibrium 
phase transition issue, including a rather accurate characterization of the as- 
sociated transition order. In the following section we set down to apply the 
Lynden-Bell argument to the present case, benchmarking the theory to numer- 
ical experiments. 




Figure 5.4: (Color online) Phase space density f(8,p) for 6 = 0.05 and e = 
0.2 (top) and 0.4 (bottom). 

Before ending this section, we report observations of phase-space structures 
resulting from the Vlasov-based simulations. Two phase-space portraits are 
enclosed in Fig. 5.4, and refer to different values of the energy e, respectively 
below (upper panel) and above (lower panel) the critical transition energy rel- 
ative to the selected (fixed) bo amount. When the system evolves towards a 
state at I ^ 0, then f(9,p) shows a large resonance. At variance, in the oppo- 
site regime, a hole-resonance dipole structure is observed, see also (Antoniazzi 
et al., 2005). This observation seems to suggest that the out-of-equilibrium 
phase transition materializes via a bifurcation of invariant structures, an obser- 
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vation that has recently been made for the HMF model (Bachelard et al., 2008). 
These findings will be further discussed in section 5.6. 



5.3 On the violent relaxation theory 

In his work on self-gravitating systems, Lynden-Bell suggested (Lynden-Bell, 
1967) that the collisionless dynamics governed by the Vlasov equation tends to 
maximize a fcrmionic entropy. The latter is obtained from the classical defi- 
nition, where the counting of the microscopic states, compatible with a given 
macroscopic configuration, results from a combinatorial calculation, and is sensi- 
tive to the underlying Vlasov dynamics. The method was successfully employed 
in the study of the HMF model (Antoniazzi et al., 2007c; Chavanis, 2006c; An- 
toniazzi et al., 2007c) and also applied to predict the quasi-stationary amplitude 
of the FEL wave (Barre et al., 2004; Curbis et al., 2007). In these works, how- 
ever, the analysis just focused on the unstable regime (7^0): no attempt was 
in fact made to reconcile it, with the high energy homogeneous state, via the 
phenomenon of out-of-equilibrium phase transitions. 

In the following we shall review the main steps of the derivation of the violent 
relaxation theory, applied to the FEL setting. Starting from a waterbag, the 
entropy to be maximized (see section 2.3.3 for details) is: 



(5.5) 



where /o is specified in (5.3) and / is the coarse-grained distribution function. 
The maximization problem is subject to constraints on normalization, total 
momentum and energy. 

Following Barre et al. (Barre et al., 2004; Barre, 2003), maximizing the 
functional (5.5), results in the following set of equations 



/o-)= J de cf q (Cx) = i, 



f ^= J d9 sin0 QF^x) = A\ 



fo W~s J deCF 2 (Cx) = e+*-A\ (5.6) 
where £ = exp (— 2A(3sm9). The functions F and F 2 are defined as follows: 

Fo(y) = r (5.7) 

J -co 1 + y e 2 

°° (5.8) 
+ y e 2 

(3 and x are (rescaled) Lagrange multipliers and ultimately stem from the 
conservation of mass, momentum and energy. A, (3 and x are calculated by 



r°° i 
F2(y) = / - 

J-co 1 
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solving Eqs. (5.6) numerically via a Newton-Raphson method (Press ct al., 
1996). The resulting (real) value of A is expected to return an estimate of the 
laser intensity at (Vlasov) saturation, I — A 2 , while f(0,p) is: 



f(0, P ) = fo 



1 



(5.9) 



\+X e P(p 2 /2+2Asine+A 2 p+A' l /2) ' 



A useful simplification occurs for A = (namely I — 0) and the optimization 
problem (5.6) reduces to: 



There exists a value of x which solves the above equation for any choice of 
AO. The homogeneous state is a stationary solution of the Lyndcn-Bcll entropy 
and thus a potentially attractive state of the Vlasov dynamics. Additional 
inhomogeneous solutions (A ^ 0, or, equivalently, 1^0) might however emerge 
from investigating the full system (5.6). The homogeneous and inhomogeneous 
solutions will be referred to as to LBO and LBA, respectively (see Fig. 5.5). 
The forthcoming discussion will focus on how to discriminate between the two, 
and eventually predict the asymptotic fate of the system. 



5.4 Interpreting the out-of-equilibrium transi- 
tion as a dynamical switch between two so- 
lutions of the theory 

Let us focus first on the LBA solution. In Fig. 5.6, we report the value of 
I as it follows from Eqs. (5.6), for different choices of the energy and initial 
bunching. The predicted intensity is shown to decrease when the energy gets 
larger but no transition is observed, in contradiction with the results of our 
numerical simulations. As previously stressed, the homogeneous LBO state is 
also solution of the optimization problem (5.6 ) and could in principle prevail over 
the former. To shed light on this issue, we calculated the entropy values Sa and 
5*0, associated to LBA and LBO, respectively. Results of the computations are 
shown in Fig. 5.7, where the dependence on the energy e is monitored for various 
choices of bo- Surprisingly, and at odds with what happens for the HMF model 
(Antoniazzi et al., 2007c), Sa is always larger than So. The two curves do not 
cross each other and the LBA configuration is entropically favored. Let us note 
that for higher values of e, I decreases and the two solutions get close to each 
other, also from the point of view of the entropy, without crossing. 

The observed transition could possibly stem from a purely dynamical mech- 
anism, and this would justify the discrepancy between the simulation output 
and the statistical prediction. More specifically, we here argue that, depending 
on the selected initial conditions, the system explores a local basin of attraction 




(5.10) 
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Figure 5.5: (Color online) Examples of f(0,p): the waterbag (top-left), the LBO 
solution (top- right), the LB A solution (bottom- left) and their corresponding 
velocity distribution function integrated over space. 



and struggles to find its way to the deputed, global maximum of the entropy. 
To clarify this point, we focus on a single numerical simulation, assuming the 
system to be initialized in a LBO state. The dynamics can progressively take 
the system towards the LBA configuration, leading to the maximization of the 
Lynden-Bell entropy. The opposite is not possible, and a simulation started in 
the LBA state will certainly not evolve to the LBO. However, dynamical effects 
might be also at play and interfere with the ideal situation here schematized, 
by virtually blocking the system in the neighborhood of an initially assigned 
LBO configuration. Is this the correct scenario? And how to explain the ob- 
served transitions that instead relate to the waterbag initial condition? These 
issues are addressed in the following, where the stability of LBO and LBA is 
investigated via direct Vlasov simulations. 

In Fig. 5.8, the dynamical evolution of the intensity I is depicted, for three 
different classes of initial conditions, relative to the same choice of e and bo. 
The LBA is indeed stable, no deviation from the initial configuration being 
observed as an effect of the Vlasov dynamics. Conversely, the LBO condition 
proves unstable, and the intensity converges towards an oscillating plateau. 
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Figure 5.6: (Color online) The saturated intensity I for the (inhomogeneous) 
LBA solution of system (5.6) plotted as function of the energy e. Different 
curves refer to distinct values of the initial bunching 6q (see legend). 



Interestingly, the LB0 and waterbag (WB) evolutions are qualitatively similar, 
and, moreover, display the same average asymptotic value for the intensity /. 
Even more important, the asymptotic value corresponds to the LBA (maximum 
entropy) solution. 

To further elucidate the analogies between LB0 and WB, and also clarify the 
issue of stability, we performed an extensive campaign of simulations, aimed at 
generalizing the results of Fig. 5.8. The results of the investigations are reported 
in Figs. 5.9 and 5.10, where I is depicted versus e, for two selected b amounts. 
The LBA initial condition is always stable under Vlasov dynamics. It is in fact 
a global maximum of the Lynden-Bell entropy, which, in this respect, proves 
adequate to describe the system at hand. The observed LB0 evolution is by far 
more complex. For large values of the energy, LB0 is stable. The stability is 
eventually lost when reducing the energy parameter: A transition materializes 
and the LB0 evolves towards the LBA state. In the vicinity of the transition, 
the LB0 initial condition approaches an asymptotic configuration, which slightly 
differs from the LBA one, and possibly results from a balance between opposing 
dynamical strengths. The WB evolution mimics that of LB0, the two curves 
returning within a pretty close correspondence. In practice, during a short tran- 
sient, an initially bunched WB expands (almost) ballistically, the particles being 
essentially transported by their own initial velocities, so visiting the whole inter- 
val [— 7r;7r]. The obtained distribution can be approximated by a homogeneous 
LB0 state (the field has not yet developed, its intensity being effectively negli- 
gible), which in turn explains the observed correspondence. We shall however 
emphasize that this mechanism applies to relatively small bo amounts. This 
fact is testified in Fig. 5.10: The agreement between LB0 and WB evolution 
is shown to worsen, when compared to that of Fig. 5.9 This is understood 
as follows: Starting from a high degree of bunching, the induced field opposes 
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Figure 5.7: The Lynden-Bell entropy calculated respectively for the LBO (plain 
line) and LBA (dash-dotted line) solution. Here bo = 0.00, 0.50, and 0.90 (resp. 
panels 1,2 and 3). 



the natural ballistic contribution, by further enhancing the tendency to form a 
coherent clump of particles. 

In summary, our calculation returns two stationary points of the Lyndon- 
Bell entropy. The first, which we termed LBA, corresponds to a inhomogcneous 
(laser on) configuration and it is a global maximum of the entropy. The second, 
labelled LBO, is homogeneous (laser off). For sufficiently large energies, the 
system can be locally trapped in the vicinity of the LBO. This happens if the 
system is initiated close enough to a LBO state, as e.g. in the case of WB 
with moderate bo values. For smaller e, the LBO loses stability and the system 
departs towards the entropically favored LBA state. 

Having detected no additional stationary points of the fermionic entropy, 
other than LBO and LBA, we interpret LBO and LBA as a saddle point and a 
global maximum, respectively. While LBO could also be in principle a global 
minimum, this hypothesis is invalidated by the fact that there exists (at least) a 
function /, compatible with the system dynamics, which yields a lower entropy 
value. Consider in fact : 

This is not a stationary point of Lyndon-Bell's entropy but represents one of 
the admissible stationary solutions of the dynamical system (5.2), in its coarse- 
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Figure 5.8: (Color online) The intensity / as a function of time, for three different 
choices of initial condition: WB (plain line), LBA (dashed line) and LBO (dotted 
line). All conditions refer to 60 = 0.05,e = 0.10 

grained perspective as implied by the theory. The Lynden-Bell entropy associ- 
ated to Eq. (5.11) can be straightforwardly computed via Eq. (5.5) and it is 
found to be always smaller than So, the value associated to the LBO configura- 
tion, for all 60 and e. Based on the above, and as previously anticipated, we can 
exclude the possibility for LBO to be a global minimum. Following our deductive 
reasoning, we hence suggest that LBO is instead a saddle-point and the observed 
transition is consequently interpreted to stem from a local modification of LBO 
stability properties or morphological characteristics (e.g. width/flatness of the 
stability basin). A detailed analytical characterization of LBO stability is at 
present missing and could eventually help clarifying the underlying scenario. 

5.5 Summary for the out-of-equilibrium phase 
transition 

Using the case of a FEL as a paradigmatic example, we described the out- 
of-equilibrium dynamics of a mean-field model for wave-particle interaction. 
Our numerical investigation moves from the Vlasov version of the model, which 
rigorously applies to the continuous limit and is believed to constitute the correct 
interpretative framework to elaborate on the QSS peculiarities. Working within 
this context, and assuming a specific class of initial conditions, we identify a 
switch between different macroscopic regimes. Such a transition is ultimately 
controlled by the nominal energy value and by the initial particles bunching, in 
qualitative agreement with what was previously observed for the HMF model 
in Ref. (Antoniazzi et al., 2007c). 
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Figure 5.9: (Color online) The saturated intensity I is plotted as function of the 
energy e, for 60 = 0. WB, LB0 and LBA initial conditions are considered, see 
legend. 
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Figure 5.10: (Color online) The saturated intensity / is plotted as function of 
the energy e, for bo = 0.50. WB, LB0 and LBA initial conditions are considered, 
see legend. 

The Lynden-Bell violent relaxation theory is developed with reference to 
the FEL setting to quantitatively substantiate our findings. We numerically 
characterized the stability of the stationary points of the fermionic entropy 
functional. A sudden change in their characteristic is related to the occurrence 
of the observed transition, supporting the adequacy of the violent-relaxation 
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theory to describe the states reached by the dynamical system. 

As a final comment, and beside stressing the unifying picture that is here 
brought forward, we emphasize that the transition here predicted can be in 
principle observed in real devices. We regard this as a suggestion for future 
experiments which could eventually result in a direct proof on the existence of 
QSS for wave-particle systems. 



5.6 Structures in phase space 

We come back in this section to the phenomenon of phase space coherent struc- 
tures observed in section 5.2. The study of phase space topology is common 
in dynamical systems, and has already been the subject of investigations in 
mean- field models: Bachelard et al (Bachelard et al., 2008) show the occur- 
rence of regular structures in the kinetic limit for the Hamiltonian Mean-Field 
; del-Castillo- Negrete (del Castillo-Negrete, 2002) finds a coherent dipole in the 
wave-particle interactions through numerical investigations. 

These structures prevent an exploration of the complete phase space for a 
given energy, preventing the statistical theory of Lynden-Bell to match the re- 
sults of simulation. This reason has been invoked in the HMF model (Antoniazzi 
et al., 2007a; Antoniazzi et al., 2007c), the FEL (see section 5.4), non-neutral 
plasmas (Levin et al., 2008b) or gravitational systems (Mineau et al., 1990; 
Yamaguchi, 2008). 

In the present section, we aim at describing these features via numerical 
Vlasov dynamics. Thereafter, we present the perspective of integrating Poincare 
sections of particles evolving under a field (A x (t), A y {t)) prescribed by a numer- 
ical integration of the Vlasov equation. 

Ultimately, we would like to obtain a reduced model for the FEL describing 
these different structures, which would allow a better theoretical understanding. 

Whereas we restricted ourselves in the present chapter to the resonant situ- 
ation (see section 5.1), we make use in this section of the detuning parameter 5 
that quantifies the deviation from this ideal behaviour. Varying 8 unveils a rich 
phenomenology in the dynamical evolution. 

The Vlasov-wave equations (5.2) become, with this addition : 

S = -p%+2(A xC o S 0-A y su i 0)^ p , 
dA x 



Oz 
dA, 



dz 

and the corresponding energy is 



= -5A y + J dddp /cos 8, 

= SA X - J d9dp / sin 6», (5.12) 



U[f](A) = Jd8dp f (8, p)(^- + 2 (A x sin 8 + A y cos 8)j- 5 {A 2 x +A 2 y ) . 

(5.13) 
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5.6.1 The resonant case 



Setting 5 — in the FEL model, the initial condition of the waterbag type is 
defined uniquely by U and bo- 

The most common occurring structure is the macroparticle, leading to a high 
value of the field. We depict it in Fig. 5.11. 




Figure 5.11: The resonant macroparticle in the FEL. U = 0.4 and bo = 0.2. 

When the macroparticle becomes less pronounced, we observe in addition a 
depression ("hole") in /, travelling in the opposite direction as the macroparti- 
cle. This situation is represented in Fig. 5.12. 

In other situations, the macroparticle is not homogeneous and presents an 
internal structure. See for instance Fig. 5.13 

These situations cannot be properly described by Lynden-Bell theory : the 
theory predicts a distribution function in phase space that depends on the en- 
ergy of the fluid element only, with a monotonously decreasing dependence (an 
illustration of the result for the theory can be found in Fig. 5.5). 

Indeed, the parameters given in Figs. 5.12 and 5.13 fall in the region of 
disagreement between the theory and the simulations. For these parameters, 
the structure present in phase space prevent the possible description of the 
state attained by the dynamics by a solution of Lynden-Bell theory. The results 
of the investigation for bo = 0.5, of application for Fig. 5.12 where U = 0.6 and 
bo = 0.5, can be found in Fig. 5.10. The evolution of the initial waterbag leads 
to a value for the intensity far below the one predicted (LBA solution) . One may 
also observe in Fig. 5.12 that a part of the initial waterbag behaves similarly to 
orbits outside of the separatrix of an equivalent pendulum and therefore spread 
uniformly on the interval [— ir;w]. This fraction of the phase space does not 
contribute to the bunching, diminishing the amplification of the wave. 
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Figure 5.12: A resonant macroparticle and a depression. U = 0.6 and b = 0.5. 




Figure 5.13: An inhomogeneous macroparticle. U — 0.6 and bo = 0.9. 



5.6.2 The detuned case 

The study performed in sections 5.1 through 5.4 is restricted by the resonant 
condition (<5 = in the model). The phenomenology of dynamical regimes is 
however richer if we relax that condition. We present, as an illustration, some 
numerical results of interest because of the intrinsic difference they show with 
respect to the resonant case. 

Figures 5.14 and 5.15 display a dipolar structure and a tripolar structure 
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Figure 5.14: A dipolar structure in the detuned case. U = 4.2 10~ 4 , bo = and 
5 = -1.40. 
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Figure 5.15: A tripolar structure in the detuned case. U = 4.2 10 4 , b = and 
S = -1.05. 



respectively. Similarly to what happens in section 5.6.1, these results cannot be 
described by a solution of Lynden-Bell's theory. 
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5.6.3 Poincare sections 



A common technique to investigate nonlinear dynamics is the Poincare section. 
It consists in observing the variables of a full dynamical system at discrete times 
in a reduced space. The times at which we collect these position are defined in 
terms of the variables of the complete system. 

We will use a criterion based on A(i) to collect data in the FEL. This 
technique has been applied recently to the HMF model, in a A-body setting, in 
order to understand the collective behavior in the thermodynamic limit. 

We present in this paragraph a technique that allows a quantitative improve- 
ment for the study of phase space. We consider a system of particles subjected 
to an external forcing A(t) = (A x (t), A y (t)) where the forcing is taken as the 
time series resulting from a Vlasov simulation. The positions and velocities of 
a set of particles are collected at each passage of / = A\ + A 2 y through its mean 
value. 

The use of Poincare sections in the study of the HMF and FEL models has 
already been used in a full A-body setting (Morita and Kaneko, 2006; Bachelard 
et al., 2008). Both papers indicate that the dependence on the number of parti- 
cles A is significant. Our approach is original and removes the dependence on A 
from the start. It allows the use of Poincare sections in a new context, making 
its power of interpretation available to understand the dynamical properties of 
the Vlasov equation. 

The Poincare section displays structures with significant details and an anal- 
ysis of the trajectories allows us to identify periodic orbits. Figures 5.16 and 
5.17 display two qualitatively different regimes, indicating a bifurcation between 
different phase space structures. 

Figure 5.16 contains an elliptic-like structure, similar to the macroparticle 
observed earlier. Figure 5.17 on the other hand shows a period-2 orbit : around 
(0,p) ~ (5,-1), a single particle finds itself successively near (4.5,-2) and 
(5.1,-0.8). 

A wide range of behaviors are found in the system while varying S, including 
orbits of higher periods. Work is in progress to extend the results towards the 
understanding of the bifurcation at the origin of the period-2 and higher period 
orbits. 
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Figure 5.16: Poincare section for the FEL. The particles are subjected to the 
force field of a Vlasov simulation, for parameters 60 = 0, U — 4.2 10 -4 and 
S = —0.5. The color code is arbitrary. Phase space displays an elliptic fixed 
point at (0,p) « (5, —1). 




Figure 5.17: Poincare section for the FEL. The particles are subjected to the 
force field of a Vlasov simulation, for parameters 60 = 0, U = 4.2 10~ 4 and 
S = —1.0. The color code is arbitrary. Phase space displays an orbit of period 
2 at (9,p) « (5,-1). 
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Chapter 6 



Relaxation in the 
Hamiltonian Mean-Field 
model 



Vlasov dynamics exhibits incomplete relaxation in a number of situations. We 
propose an analogy with the mixing in fluids and apply it to the study of Vlasov 
dynamics in phase space. We find that stretching and folding structures appear in the 
evolution of the Hamiltonian Mean-Field model. In order to quantify the development 
of these structures, we compute the contour of the fluid in phase space through high 
accuracy simulations. 



The statistical mechanics of systems with long-range interactions is known 
for simple models with good detail (see for instance (Campa et al., 2009)). The 
approach to equilibrium remains however a challenging aspect in these systems. 

For instance, given appropriate initial conditions, a system may not reach 
equilibrium in the thermodynamic limit. It is instead trapped in a so-called 
quasi-stationary state (QSS) whose lifetime grows with the number of particles 
(Antoni and Ruffo, 1995; Latora et al., 1998). 

A most remarkable achievement is the prediction of the state of the system 
in these QSSs thanks to Lynden-Bell's theory (Lynden-Bell, 1967; Antoniazzi 
et al., 2007c). Lynden-Bell's theory bases itself on the microscopic properties 
of Vlasov dynamics and on the assumption that there is a mechanism leading 
the system to its most probable state. The framework of the Vlasov equation 
leads to a better understanding of the dynamical properties of systems with 
long-range interactions (Yamaguchi et al., 2004; Barre et al., 2006; Jain et al., 
2007; Antoniazzi et al., 2007a). 

It is argued (Antoniazzi et al., 2007a; Levin et al., 2008b; Levin et al., 
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2008a) that the discrepancies of Lyndon-Bell's theory to successfully predict 
the outcome of a system for some ranges of initial conditions are due to purely 
dynamical effects. 

In order to detail the relaxation process taking place in phase space, we 
propose an analogy with fluid mixing. Ottino (Ottino, 1989) gives a detailed 
introduction to the dynamical theory underlying mixing in fluids. 

Taking the Hamiltonian Mean-Field (HMF) model as an example, we com- 
pute properties of the dynamics from the phase space fluid point of view. We 
identify regions in which no dynamical evolution takes place, whereas the con- 
tour of the distribution function experiences substantial deformations in phase 
space. 

We would like to add that the results presented in this chapter would greatly 
benefit to the understanding of relaxation in the gravitational sheet model. Ex- 
tensive studies of this model exist in the literature (Luwel and Severne, 1985; 
Funato et al., 1992; Yamaguchi, 2008) and the approach that we present would 
complete adequately the dynamical considerations, especially those of Ref. (Fu- 
nato et al., 1992) in which a separation of the initial waterbag into distinct 
regions is performed in phase space. 

This chapter is structured as follows: section 6.1 defines the HMF model 
and the initial conditions that we consider, section 6.2 introduces the analogy 
between fluid dynamics and phase space dynamics, section 6.3 details results 
for the simplified case of the pendulum and section 6.4 presents the results for 
the HMF model. 



6.1 Dynamical evolution of the HMF model 

The Vlasov equation describing the HMF is: 

?L + M dv if]df _ 
dt p de do dp 

V[f}(6) = 1 - M x [/] cos 9 — My [/] sin 9 (6.1) 
where M = (M x ,M y ) is defined by: 

M x [f] = J dOdp /cos0 ; M x [f} = J d9dp f sin 9 , (6.2) 

and the magnetization M is the modulus of the vector (M). Equation 6.1 is 
solved numerically with the code presented in chapter 4. 

The investigations in this chapter are performed for a waterbag initial con- 
dition defined by: 

' fo if \P\ < Ap, 
f(9,p) = { \9\< AO, (6.3) 

otherwise, 
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where AO and Ap characterize the extent of the waterbag in both directions of 
phase space. Normalization imposes /q = 4A g Ap ■ We can also characterize the 

initial condition uniquely by M = and U = + ^(1 — Mq). 

We use throughout this chapter the waterbag initial condition with U = 0.6 
and Mq = 0.2. Figure 6.1 displays the magnetization M as a function of time. 
The parameters used in Fig. 6.1 are: Ng = N p = 512 and At — 0.1. M(t) 

0.45i 1 , , 1 , , , 1 




t 

Figure 6.1: The magnetization M in function of time for a waterbag initial 
condition with U — 0.6 and Mq = 0.2. 

begins with strong oscillations. After t s» 100, the oscillations become regular 
and damp slowly. 

6.2 Mixing in fluids and in phase space 

The region of phase space defined by the initial waterbag can be regarded as 
a fluid. The Vlasov equation presents the form of a complicated advection 
equation deforming the said fluid. This deformation leads to a filamentary 
structure even in the simplest setup of free streaming (see section 2.2.5). 

The Vlasov equation possesses the property of keeping levels constant: the 
one-particle probability distribution function, starting with and /o as the only 
two possible values, will keep and /o as the only accepted values during the 
whole time evolution. The conservation of the normalization in the course of 
time implies that the surface occupied by the fluid in phase space also remains 
a constant. The perimeter of that region can grow in time as the shape of the 
fluid evolves and develops filaments. We expect it to behave similarly to a 2D 
fluid submitted to an external perturbation, a situation giving rise to mixing in 
fluid dynamics (Ottino, 1989). 
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The length of the contour is monitored in the course of time and serves as a 
quantitative indication of the deformation of the waterbag. 

For practical purposes, we define the perimeter Pf(t) as the length of the 
contour of f(0,p) whose height is fo/2. This allows a convenient numerical com- 
putation in simulations through the use of the CONREC Fortran subroutine by 
Bourke (Bourke, 1987). We implemented this routine in the Vlasov simulation 
code presented in chapter 4 

6.2.1 The free-streaming case 

An analytical computation of Pf(t) is possible in the case of free streaming. 
Pf(t) grows linearly in time, except for a small deviation for short times: 

P f (t) = 4A6» + 4Ap V / l + t 2 . (6.4) 

This can be readily checked through numerical simulation. Figure 6.2 dis- 
plays a perfect agreement between Eq. (6.4) and the numerical computation. 




Figure 6.2: Evolution in the free streaming case of the perimeter Pf(t) for an 
initial waterbag with AO = 2 and Ap = 1. The simulation matches exactly the 
expected theoretical behaviour. 

We expect that a linear behavior of Pf (t) will characterize a regular dynam- 
ics, as in the case of integrable systems. 

6.2.2 Dynamics in the presence of interactions 

The contour of the waterbag experiences complicated deformations in phase 
space. We give a sketch of this idea in Figs. 6.3 and 6.4 
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Figure 6.3: Illustration of the deformation of a waterbag. Between the two 
arrows, the contour experiences stretching, giving rise to a growth of Pf(t). 




Figure 6.4: Illustration of the deformation of a waterbag. After the stretching, 
the contour folds on itself. 



The repeated action of stretching and folding modifies considerably the orig- 
inal waterbag. It is recognized as a signature of a chaotic flow in fluid dynamics 
(Ottino, 1989). 

In Fig. 6.3 , the deformation pointed by the arrows leads to a stretching of 
the contour. The distance between two nearby points on the contour increases, 
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causing the length of the contour to increase. If the separation between these 
points grows exponentially, the behavior of Pf(t) will contain an exponential 
contribution : 

P f (t) cx e xt . (6.5) 

If the average magnetization becomes a constant after some time tc, we expect 
that the perimeter acquires a linear behavior after tc, similar to the case of the 
pendulum (see below and in Fig. 6.5). 

Pf(t) being accessible in simulations, we dispose of a method to quantify 
the deformation of the waterbag in the course of time. The observation of 
an exponential behavior would then provide the a posteriori validation of our 
hypothesis. 

6.3 The pendulum 

The Hamiltonian for the HMF model is very similar to the one of a pendulum. 
In the pendulum, the field is set to a constant value while in the HMF model 
it depends on the positions of the particles. This similarity is exploited in 
chapter 7 with great detail and we discuss here only what is needed to allow the 
interpretation of the fluid analogy. 

For the initial waterbag given by U — 0.2 and M = 0.6, M y = at all times. 
The dynamics of the HMF model is thus similar to the one of a pendulum whose 
field intensity varies in time, given by the following Hamiltonian: 

ff = Ef- M wE cos ^ • ( 6 - 6 ) 

i i 

The phase space of a pendulum is divided by the separatrix in two qualita- 
tively different parts : particles are trapped if \p\ < otherwise they travel 
periodically with a velocity of constant sign. 

The case of a constant field intensity M(t) — M(0) already induces a de- 
formation of the initial waterbag. However, the separation between two nearby 
points of the contour remains bounded linearly. Figure 6.5 displays the evolu- 
tion of the perimeter in the case of the pendulum alongside with a linear fit. 



6.4 Structures and relaxation in phase space 

We now turn to the study of the dynamical evolution of the Hamiltonian Mean- 
Field (HMF) model. 

The evolution of the initial waterbag is displayed in Fig. 6.6. The central 
part of the fluid remains constant (the white region), experiencing a rotation in 
phase space. The deformation appears only on the contour, with the formation 
of filaments. The length of the filaments increases in time while their thickness 
decreases. 
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Figure 6.5: Evolution of Pf(t) in the case of a pendulum with a constant field 
intensity. "2k" is a Vlasov simulation with parameters Ng = N p — 2048 and 
At = 0.1, the initial condition is a waterbag with AO w 2.6 and Ap « 0.85. 
Pf(t) is approximated by a linear fit ("linear" in the legend) around which 
there are small variations. 



In the region containing the filaments, a coarse-grained probability distribu- 
tion function f(6,p) can take any value between and fo- The phase space is 
thus separated qualitatively in three regions: 

1. The central region where f(0,p) is constant at the value fo. 

2. An intermediate region where filaments are formed and where an effective 
mixing takes place. 

3. The remainder of phase space, where f(9,p) is equal to 0. 

The computation of Pf (t) is performed in the course of the simulation. The 
result is displayed in Fig. 6.7 and exhibits an exponential behavior (the scale 
of the figure is logarithmic). A fit ("exp" in the legend) to an exponential is 
superimposed on the figure. After a short time, the level of detail required to 
follow accurately the filamentary structure exceeds the accuracy of the numerical 
grid. The exponential behavior is then lost. 

Successive simulations with higher accuracies have been run and are also 
displayed in Fig. 6.7. The time during which Pf(t) agrees with the exponential 
fit increases slightly with the grid accuracy. We expect successive increments 
of the number of grid points to lengthen that time of validity. This could be 
checked with either a parallelization of the simulation code or with the help of 
alternative algorithms. 
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Figure 6.6: Evolution of the waterbag in phase space, in the case of the HMF 
model. Parameters of the simulation: Ng = N p = 512, At = 0.1, U = 0.2 and 
Mo = 0.6. The usual pseudo-color representation is completed by the contour 
of f(9,p) for the value fo/2. The central part of the waterbag remains constant 
but the contour of the waterbag displays filamentation. 

We display for the highest number of grid points a zoom on the contour lines 
in Fig. 6.8. A folding structure is clearly apparent in the upper region of the 
figure. This is a feature that is also found in fluid dynamics (Ottino, 1989). 

A study of the iV-body dynamics of the HMF model (Sguanci et al., 2005) 
reveals that the dynamical evolution for high value of N (10 6 ) does not lead to 
fractal structures. The reference (Sguanci et al., 2005) uses a singular initial 
condition that is not tractable by numerical simulation for the Vlasov equation, 
preventing a direct comparison of numerical results. 

6.5 Summary 

We provided in this chapter original contributions to the understanding of the 
approach to equilibrium in collisionless systems. In the Vlasov equation for the 
Hamiltonian Mean-Field model, we considered an analogy with fluid dynamics. 
The analysis of phase space shows a separation in different regions, one of which 
displaying mixing. 

The presence of stretching and folding structures in phase space has been 
observed in numerical simulations. The deformation of the contour of an initial 
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Figure 6.7: Evolution of Pf(t) for the HMF model, same conditions as in Fig. 6.6 
but with different numerical accuracies. The number of grid points in runs "2k" , 
"4k" and "8k" is Ng = N p = 2048,4096 and 8192 respectively, "exp" is an 
exponential fit of the region < t < 12 

waterbag has been quantitatively monitored in the simulations and supports 
the hypothesis of an exponential elongation linked to the stretching and folding 
phenomenon. 

A pursuit of this work is possible through the use of more advanced compu- 
tational techniques. In addition, a systematic study of the dynamics is necessary 
in order to perform useful comparisons to the existing body of literature on the 
relaxation of collisionless systems. 
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Figure 6.8: A zoom on a region of phase space for the HMF model. Parameters 
of the simulation: Ng = N p = 8192, At = 0.1, U = 0.2 and M = 0.6 (same 
conditions as in Fig. 6.7). The value of f(9,p) in the shaded region is /o and in 
the white region 0. The blue lines indicate the contour computed by the CON- 
REC subroutine. The upper part of the figure displays a series of intertwined 
contour lines resulting of the stretching and folding mechanism. 
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Chapter 7 



Dynamics of the 
Hamiltonian Mean-Field 
model on the basis of the 
equivalent non-interacting 

system 



With the aim of understanding the mechanisms of violent relaxation and of out-of- 
equilibrium phase transitions, we introduce a Vlasov equation for a set of uncoupled 
pendula and the corresponding Lynden-Bell theory. The asymptotic evolution of the 
system is found to match the prediction of the theory for a range of parameters. We go 
on by computing under an ergodic-like hypothesis the exact solution for the asymptotic 
evolution of the system. The stationary regime of this system is identical to the one of 
the Hamiltonian Mean-Field model, allowing us to improve our understanding of this 
latter model. We find an out-of-cquilibrium phase transition under a self-consistent 
constraint. 



The evolution of the Vlasov equation displays complex nonlinear dynamics 
in the presence of interactions. We propose to study a Vlasov equation that does 
not include interactions, providing an integrable example of Vlasov dynamics. 

The steady state of the Hamiltonian Mean-Field (HMF) model corresponds 
to a constant value of the magnetization vector. The equations of motion are 
then equal to the ones of the set of pendula. This similarity motivates the inves- 
tigation of the Vlasov dynamics for a set of uncoupled pendula. Lynden-Bell's 
theory aims at predicting the outcome of collisionless evolution and supposes 
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the existence of a mixing mechanism in the microscopic dynamics. We develop 
it for our intcgrable system to test that particular hypothesis. 

The investigation is pushed a step further by the computation of an exact 
solution to the asymptotic evolution of the set of pendula. Imposing a self- 
consistency condition explicitly specializes our study to situations comparable 
to the HMF model. The idea to study the dynamics of a set of uncoupled 
pendula has been put forward recently in the literature (Leoncini et al., 2009; 
Firpo, 2009) to construct a class of stationary solutions for the HMF model. 
Our approach differs in the sense that we detail the computation explicitly with 
the aim to understand the occurrence of a phase transition. 

This chapter is structured as follows: the system of pendula and its asso- 
ciated Vlasov equation are introduced in section 7.1. Lynden-Bell's theory is 
developed in section 7.2 and the exact solution to the asymptotic dynamics is 
given in section 7.3 



7.1 Vlasov equation for the pendula 

The system we consider consists of uncoupled pendula subjected to a fixed 
external field Ti : 

N o N 

i=i j=i 

where the pj are the momenta, the 9j the positions and Ti is the external field 
applied to the system. We introduce the bunching parameter 



= (b x ,b y ) = ^ ^cos^^sin^ (7.2) 




and its modulus b — ^jb% + b y . In this system, a symmetric initial condition 

with by = leads to b y — at all times. We consider only such situations 
in this chapter. The interaction part of the Hamiltonian, rescaled by N, can 
be rewritten as = —Tib. We also introduce the one-particle Hamiltonian 

Hi{6, P ) = £-Ticos6. 

The bunching parameter b is computed with the same expression as the 
magnetization in the HMF model but does not enter in the force acting on each 
particle. This fundamental difference implies that the dynamics of system (7.1) 
is integrable. The value of b indicates the homogeneity of the system: 6 = 
implies that the system is homogeneous and 6^0 indicates an inhomogeneous 
system. 6=1 means that all particles are located at 6 = 0. We make use of 6 
to quantify the state of the system or track its evolution in time. 

Stationary regimes of the HMF model (constant value of the magnetization 
M) possess the same equations of motion as the ones of system (7.1) : 

{ p] =-Wsin^ (? ' 3) 
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allowing a comparison of QSSs of the HMF model with the ones of the uncoupled 
system. 

In the limit N — > oo, we describe the evolution of the system from the kinetic 
point of view, using the one-particle probability distribution function (PDF) in 
phase space. The PDF, f(6,p;t), gives the probability f(6,p;t)d6 dp to find at 
time t a particle in the phase space region centered on (6,p) and of measure 
d6 x dp. We write the following Vlasov equation 



0/ Of dV{6) df _ 

dt p de de dp 



(7.4) 



where V(9) = —ficosO is the potential. This formulation, although similar 
to other Vlasov equations, is different in the sense that it does not include 
interactions, but only an external field. 

A consequence of the fact that the potential does not depend of / is the time 
independence of the energy level distribution function defined as: 



p(e) 



-! 



MdpffrpWe-HMp)) 



(7.5) 



Two qualitatively separated regimes take place whether H is equal to or 
not. The evolution of b(t) for these two situations is given in Fig. 7.1 and we 
depict the evolution of the waterbag initial condition by the Vlasov dynamics 
in Fig. 7.2. 
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Figure 7.1: Evolution of the bunching parameter b in function of time for the 
two situations % — and H^0. 

We notice from Figs. 7.1 and 7.2 that in the free-streaming case (H = 0), 
the bunching parameter b goes to zero, corresponding to a spread of the initial 
waterbag to the whole [— 7r; ir] interval, whereas in the T-L ^ case a finite value 
of b is attained. The spreading of the initial waterbag still occurs but follows 
equal energy lines in the phase space of the pendulum instead of straight lines. 
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Figure 7.2: Evolution of a waterbag initial condition in phase space. The left 
column represents the H — case, from top to bottom chronologically. The 
right column represents the % ^ case. 



The complete range of waterbag initial conditions are given equivalently by 

f-Hb and b = s -^ 



the parameters AO and Ap or U = Hbo and bo — 81 . The waterbag is 



defined as 

r /o if h < a p , 

f(0,p) = \ \0\< A6, (7.6) 

[ otherwise. 

with f = (4A6>A P r\ 

The H 7^ regime allows one to rescale the system in order to explore the 
possible regimes more easily. We perform the following changes: 

H -> H' = H/H, (7.7) 
Pi^ Pi=Pi/VU, (7.8) 
Bi -> fl{ = ^ (7.9) 

and obtain a system in which % is effectively set to 1. 

The rescaled system will be used when appropriate with explicit mention. 
For the comparison with HMF stationary states, we however must find a solution 
where b — Ji, which is trivial under the inverted rescaling. As 9 is not affected 
by the rescaling, the bunching is the same in both systems and to find % = b, 
we only have to find a set of AO and Ap for which the final b is the desired 
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one and rescale the system to adjust T-L. We remark however that homogeneous 
situations with % = are not amenable to the rescaling. 



7.2 Lynden-Bell's theory for the pendula 

Lyndcn-Bcll devised a theory whose aim is the prediction of the outcome of the 
dynamics of the Vlasov equation ( Lynden-Bell, 1967). Its original purpose was 
to explain observations made on the light distribution from elliptical galaxies. 

The theory is based on the maximization of an entropic functional under ap- 
propriate macroscopic constraints (see section 2.3). The underlying hypothesis 
is that the dynamical evolution of a collisionless system possesses a mechanism 
to explore the ensemble of states compatible with the constraints. 

The set of uncoupled pendula allows, in principle, a formulation of the theory. 
Its formulation is indeed very similar to the one for the Hamiltonian Mean-Field 
(HMF) model as it describes only stationary regimes. The conservation of the 
energy distribution function implies that relaxation is not likely to occur in the 
system of pendula. We propose to check the possible application of the theory 
to this system in order to deepen our understanding of collisionless dynamics. 



7.2.1 Formulation of the theory 

We write here the expression of the entropy as applicable to a two-step initial 
one-particle probability distribution function (1-PDF) f(0,p) of heights and 
/o: 

s (f) = - J dp de 1 1 ' J ■ f < 1 ^ ' ' 1 



fa fo 



fo 



In 1 



fo 



(7.10) 



where f(0,p) is the coarse grained 1-PDF. This coarse grained distribution can 
take values in the whole interval [0 ; fo]. The solution / of the maximization 
problem for the entropy (7.10) is a prediction for the asymptotic outcome of 
the dynamics. It is a prediction of the 1-PDF if one neglects the microscopic 
filamentary structures that develop in phase space. 

Under the constraints of mass, momentum and energy conservations, the 
optimization of the entropy (7.10) for an initial waterbag, 



f(6,p) = f if \e\ < A6 and \p\ < Ap 
= , else 



(7.11) 



leads to the following f(6,p): 



fo 



4 A6 Ap 



(7.12) 



/(M = /o 



e -/3(p 2 /2-«cos0)-^ 
1 _|_ e -/3(p 2 /2-«cos0)-Ai 



(7.13) 
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where /i and ft are Lagrange multipliers associated respectively to mass and 
energy conservations. The resulting set of equations to be solved is: 



2ff 



^§Jd9 e^ cose F (xe^ cosfl ) = 1 
$fefd0 e^ cos9 F 2 (xef mcos9 )+H(l-b) = U 
£}§jd9 cos0 e^ cos0 F o (a;e^ cose ) = b 



(7.14) 



where x = e b is the previously defined bunching parameter and the F n 's are 
defined as: 



We solve the system of equations (7.14) with the Newton- Raphson method 
(Press et al., 1996). 

The values of b and T-L are not determined a priori. We solve the system of 
equations under the additional constraint that 



that we call the self-consistency relation. Taking into account Eq. (7.16), we set 
A8, Ap and the numerical procedures gives b — Ti , f3 and x. 

We remark that the system of equations ( 7.14) is similar to the equivalent 
problem for the HMF (Antoniazzi et al., 2007b). The value of U needs to be 
adapted to accommodate the definitions of both systems. 



A first comparison is drawn for initially bunched systems. AO is set to 1.66 and 
Ap is varied. Eqs. (7.14) are solved with the self-consistency relation b = %. A 
simulation of the Vlasov equation for the pendula with the given A9, Ap and 
the computed T-L provides a test of the validity of the theory. 

The stirring causes variations in b that diminish quickly as the initial wa- 
terbag evolves in phase space. This behaviour is displayed in Fig. 7.3. The value 
attained by b is in all situations in perfect agreement with the predicted one. 
This first result on the simulation of the pendulum bears an important significa- 
tion with respect to Vlasov dynamics: the stirring mechanism by itself plays an 
important role in the evolution of the initial condition towards a "most probable 
state" as expected by Lynden-Bell's theory for a system without collisions. 

Fig. 7.4 summarizes this result for a wider range of Ap's. The points in- 
dicating the result of Vlasov simulations match everywhere the theory. This 
is trivial when Ti is equal to 0: the system evolves under free streaming only 
and goes to a homogeneous state. It is more unexpected when ?{/0 because 
Lynden-Bell's theory, as statistical mechanics, is based on the assumption that 




(7.15) 



b = H 



(7.16) 



7.2.2 A6 = 1.66 
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Figure 7.3: For AO = 1.66, evolution of the bunching parameter in Vlasov sim- 
ulations in the course of time for AO — 1.66 and Ap = 0.41, 0.72, 1.03 and 1.34 
for panels (a),(b),(c) and (d) respectively. The corresponding value predicted 
by Lynden-Bell's theory is plotted for reference. 



the complicated dynamical evolution of the system allows an effective explo- 
ration of the states allowed in its state space. This is obviously not the case 
when the dynamics is integrablc. 

An additional result visible in Fig. 7.4 is that for growing Ap a decrease of 
the resulting b is observed. It is followed by an abrupt transition when the b ^ 
solution ceases to exist. Below the transition, both solutions to Lynden-Bell's 
theory exist, and the choice between any of the two can be made. 

To understand the effectiveness of our prediction, we detail the results of 
Fig. 7.4 with a comparison of the marginal distributions. Fig. 7.5 displays p(6) 
for some of the simulations. Lynden-Bell's theory clearly does not take into ac- 
count the filaments, which leads to a disagreement strongly visible only in panel 
(a). Else, the predicted p matches very closely the result of the simulations. 

A complementary approach is given by the velocity marginal <p(p) given 
in Fig. 7.6. There is still a marked effect of the filamentation in panel (a) 
(Ap = 0.41) but the most striking difference is in panel (d) (Ap = 1.34). The 
predicted field is and the velocity distribution is left unchanged from the 
original waterbag. The discontinuous character of the waterbag is completely 
ignored by Lynden-Bell's theory. In the system of pendula, p(e) does not change, 

2 

which translates for H = in the invariance of <p(p) (because e = 

The observed agreement is best in panel (c) (Ap = 0.72). In that situation, 
we note that the final value of b is very close to the initial one. While we find 
that the theory predicts with a good accuracy the outcome of the simulations, 
it makes sense that it works optimally when no strong changes are required 
between the start of the simulation and its end. More detail on this particular 
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Figure 7.4: For AO = 1.66, evolution ol the bunching parameter predicted by 
Lynden-Bell's theory when b — % in function of Ap. Data points correspond to 
Vlasov simulations with identical Ad, Ap and H. 

aspect will be given in the next section, where a stronger disagreement between 
the theory and the simulation motivates an deeper analysis. 

7.2.3 A6 = tt 

We now turn to homogeneous initial conditions. They differ from bunched ones 
in the sense that, in order to attain a finite bunching value, the system goes from 
a homogeneous state towards a qualitatively different inhomogeneous state. 

Figure 7.7 shows the dependence on Ap of the self-consistent bunching for 
A6 = 7r, for Lynden-Bell's theory and Vlasov simulations. For sufficiently high 
Ap the bunching is zero and the system exactly evolves under free streaming. 
For lower values of Ap, the bunching is finite but the two approaches do show 
a quantitative agreement. 

At variance with what happens for A9 = 1.66, there is no satisfactory agree- 
ment anymore between the results of Lynden-Bell's theory and the simulation. 
This appears in the panel (a) of Figs. 7.8 and 7.9 p(9) and ip(p) show a disagree- 
ment between the theory and the simulation, corresponding to the left region 
in Fig. 7.7 (low values of Ap). A slight disagreement is also observed in panel 
(a) of Fig. 7.5 , but the resulting value of b is similar to the one predicted by the 
theory in that situation. 

A quantity that we have not explored yet is the energy probability distri- 
bution function (energy PDF) p(e). For the system of uncoupled pendula, p(e) 
does not vary in time, which should prevent relaxation from happening. 

We display in Fig. 7.10 the distribution p{e) : directly computed on the 
simulation grid and compared with the distribution predicted by Lynden-Bell's 
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Figure 7.6: The p marginal, <p, for AO = 1.66 and Ap = 0.41, 0.72, 1.03 and 1.34 
for panels a,b,c and d respectively. The filamentary character is strong in panel 
a. Panel d displays a significant deviation from the theory for ip whereas it was 
exact for p (see Fig. 7.5). 101 
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Figure 7.7: For A6* = it, evolution of the bunching predicted by Lynden-Bcll's 
theory when b = % in function of Ap. Data points correspond to Vlasov simu- 
lations run with identical AO, Ap and %. 



The observation of p(e) sheds light on the effectiveness of Lynden-Bell's 
theory as discussed in section 7.2.3 and to a smaller extent in this section (panel 
(b) of Figs. 7.8 and 7.9). The theory and the simulation agree when p(e) is close 
enough to the energy PDF predicted by Lynden-Bell's theory. Then, stirring 
in phase space induces a spread of the initial waterbag along equal-energy lines 
in phase space, giving oscillations of b(t) followed by a relaxation to the value 
predicted by the theory (see Fig. 7.3). 

As the energy PDF is constant, we conclude that Lynden-Bell's theory is 
effective in the case of the pendula only when the initial condition is close 
enough to the energy distribution of the theory. When that condition is met, 
no relaxation is needed in order to attain the predicted bunching but only a 
repartition of the initial waterbag on equal energy lines in phase space. 

7.2.4 The role of the conserved quantities 

The constraints on the macroscopic quantities play an important role in Lynden- 
Bell's theory. The system of uncoupled pendula offers the perspective to modify 
the theory : the energy distribution function (DF) is conserved. The con- 
sequence is that any moment of the energy DF, E l [f], is expected to be a 
constant. 



Equation (7.17) provides the theory with an arbitrary number of conserved 
quantities that can be taken into account. This task is not possible in the case 



theory. 




(7.17) 
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Figure 7.8: The 9 marginal, p, for A9 = tt and Ap = 0.41,0.72, 1.03 and 1.34 
for panels (a),(b),(c) and (d) respectively. Apart from the remaining in panel 
(a), the agreement is very good. 
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Figure 7.9: The p marginal, ip, for A9 = tt and Ap = 0.41,0.72, 1.03 and 1.34 
for panels (a),(b),(c) and (d) respectively. The filamentary character is strong 
in panel (a) . Panel (d) displays a small significant deviation from the theory for 
if whereas it was exact for p (see Fig. 1708). 




Figure 7.10: The energy probability distribution function for A9 = tt (top) and 
Ap = 0.41 (bottom) and 0.72 (corresponding to panels (a) and (b) of Fig. 7.5 
respectively), computed by Lynden-Bell's theory and directly on the simulation 
grid. Both results are in qualitative agreement for Ap = 0.72 (lower panel) but 
differ strongly for Ap — 0.41 (upper panel). 

of the HMF model, where there is no a priori conserved quantity that remains 
to be included. The model that is introduced in this chapter thus presents an 
interest for the understanding of the role of the conserved quantities in Lynden- 
Bell's theory. 

7.3 Exact solution 

On the basis of the uncoupled character of the system described by Eq. (7.4), we 
propose an exact solution to the asymptotic dynamics. This solution is based 
on an ergodic-like hypothesis that is justified by phase space stirring instead of 
chaotic dynamics. 

The procedure allows to compute the bunching resulting from an arbitrary 
initial waterbag. A self-consistent requirement allows to discuss the similarity 
with stationary solutions of the Hamiltonian Mean-Field model. 

We explain the principle of the procedure in section 7.3.1 and its explicit 
computations in section 7.3.2. 

7.3.1 The ergodic-like hypothesis 

While ergodicity is usually considered to hold when a system experiences chaotic 
or mixing dynamics, we use this hypothesis in the following sense: every particle 
or fluid element experiences a motion that is constrained on its initial energy 
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manifold. The non-isochronic character of the pendulum gives different periods 
to two nearby orbits, leading to stirring of the phase space. Especially, the 
system does not come back close to its initial condition. 

The energy distribution p(e) is constant and we make the hypothesis that 
the asymptotic state of the dynamics f(9,p) is defined only on the basis of p(e). 

This is justified in our model because of phase space "stirring", that is, 
the evolution of the initial waterbag with no interaction leads the evolution 
towards a steady state. It is true on a coarse-grained level, similarly to Lynden- 
Bcll's theory. Everywhere in phase space, a cell of volume d9 x dp will contains 
after a time long enough, a filamentary structure. The value of / is still two- 
step valued (at the microscopic level) because the original waterbag is two-step 
valued. However, the volume is occupied by filaments and the average value of 
/ in the cell will be proportional to p(e). 

7.3.2 Explicit computation of the solution 

For the waterbag given by Eq. (7.6), the energy distribution is given by 

Pe(£,Pa) = 4K9Ap~ J d6d P 6 (lp 2 - e - ncose ) ■ ( 7 - 18 ) 
Carrying out the integral over p, one gets 

Pe(^Po) = 1A 2 nA [ d9 1 =. (7.19) 

ey ,m > AAOApJ ^2(e + ffcos0) V ; 

The integration over needs to be done on a domain enclosed in the original 
waterbag, or explicitly: 

< e + "Hcos0 < ^Ap 2 , \6\ < A9 and \p\<Ap (7.20) 

These conditions lead to two solutions, depending on the ordering of e with 
respect to H. For e>7i, the distribution function is given by : 

P e (e,A6,A P ) = -X- / d9 - = _ v / dB - =, 

AA9Ap J_ Ag ^{e + Ucosff) ±A0Ap J_ dl ^( e + Hcos6) 

(7.21) 

while for e < % it is 

P e (e, AO, Ap) = f° d9 - 1 -TTZ7T- dd ; 1 • 

AAOAp J_ 8a ^( e + ncos0 ) 4A8A P J_ dl ^( e + Hcos6) 

(7.22) 

Here Q\ satisfies 

e + H cos 0i = Ap 2 /2 , for e > Ap 2 /2 - H 

(7.23) 

0i = , else. 
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Similarly, 8q satisfies 



e + ft cos O = , for - ft < e < ft cos A6 

d = , for e < -ft ( 7 - 24 ) 
6<o = A6< , else. 



The graphical computations corresponding to Eqs. (7.21) through (7.24) are 
represented in Fig. 7.11 P t corresponds to the length of the full lines, whose 
boundaries 0\ and 0q are also displayed in the figure. 
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Figure 7.11: Graphical computation of Eqs. (7.21) through (7.24). The full lines 
represent the domain on which the integral is performed in order to compute 
P e for values of the energy e = 0,ft and 2ft. The boundaries 8q and Q\ are 
displayed for each energy. 

In the steady state, the (8,p) distribution is such that for any given energy, 
all the micro-states corresponding to that energy are equally probable, with the 
boundaries given by the initial waterbag no longer taken into account . Thus 
the steady state distribution P(6,p) may be expressed as 



P(0,P) 



1 P e (e(g, P ),Ag,Ap) 
4A0Ap Q e (e(0,p),A0, Ap) 



P e (e(0,p),A0Ap) , 



(7.25) 



where Q e (e(6,p), A0, Ap) is given by P e (e(0,p), AO, Ap), as evaluated by Eq. (7.22), 
except that one takes 

01 = (7.26) 
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and 



tt ife>U, 

(7.27) 

arccos(A#) else, 

in these integrals. The graphical computation of Q e is depicted in Fig. 7.11. Q e 
correspond to the length of the dotted lines. 

The value of P e may be interpreted, with reference to Fig. 7.11, as the ratio 
of the length of the full line by the length of the dotted line, computed for each 
value of e. One can observe that for small values of the energy e, the line is 
full on the whole length, therefore P e — (4A6>Ap) . On the other hand, for e 
"above" the waterbag, P e = 0. 

The 9 marginal may be evaluated by integration over p. 

Pe(6) = J dp P(9,p) = Jde P(M J = J de P e (e)^ . (7.28) 
Hence, by expressing the derivative in terms of e and 8, one finds 

/oo -, 
de P e (e) . (7.29) 

-ffcosfl V( e + H COSV) 

Similarly, the p marginal is given by 

P p {p) = Jde P e ( £ )^ . (7.30) 
Expressing the derivative in terms of e and p one finally obtains 

fP 2 +H i 

P P (P) = / de P t (e) (7.31) 

The elliptic integrals are performed numerically using the specialized routine 
of Ref. (Press et al., 1996). 



7.3.3 Check against a simulation 

The only hypothesis made in the derivation of the exact solution is that the 
filaments will be smoothed out by a simulation. We perform a Vlasov simulation 
to check the adequacy of our theory. 

The bunching shows in Fig. 7.12 strong oscillations in a first stage, followed 
by a relaxation to the predicted value. 

Figures 7.13 and 7.14 depict p and <p respectively. As expected, apart from 
the signs of filamentation, the agreement is excellent. 
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Figure 7.12: The bunching b in the course of time for a Vlasov simulation with 
initial conditions A9 = 2.0, Ap = 0.78 and H = 0.608. After a transitory 
oscillations, b relaxes to the predicted value. 



0.40 




Figure 7.13: p for the same simulation Figure 7.14: ip for the same simulation 
as Fig. 7.12. as Fig. 7.12 

7.3.4 Back to the HMF model 

In the stationary regime for the HMF model, a comparison is possible with the 
set of pendula. In order to proceed, one should evaluate the bunching 

/ + TT 
d£ cos 0P e (9) , (7.32) 
-7T 

and set b = H, obtaining a self-consistent relation. Clearly % = b = is a 
solution of this equation. However one can in principle find another solution of 
this equation with non-zero magnetization m. To obtain these results, we use 
an iterative root finding routine (Brent's method, see Ref. (Press et al., 1996)), 
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setting AO and Ap. The routine then seeks a root of f_* d6 cos 9Pg(9) — Ti, 
where Ti is the variable. 

The properties of both solutions are then cast in the HMF setting (where 
Ti = b is then called the magnetization M), allowing to use our exact results. 
Especially, the behavior of Eq. (7.32) can possibly display a transition between 
regions of parameters with two solutions (b = and b =/= 0) and regions where 
only b = is possible. Figure 7.15 represents the computation of Eq. (7.32) for 




Figure 7.15: Graphical explanation of Eq. (7.32) in the case A8 = 3.0. The 
intersection between the diagonal and the computed b gives the solution(s) to 
Eq. (7.32). For low Ap there can be two solutions (Ti = and T-L ^ 0) while 
above a threshold only one of the solutions exists. 

A6 = 3.0 and several values of Ap. The intersection between b and Ti gives the 
solution to the self-consistent condition. 

The value of Ti given by Eq. (7.32) allows the computation of the kinetic 
and interaction contributions to the energy in the initial state: 

Ap 2 Ap 2 

U=-%--W> , K = -%-, 7.33 
6 6 

and in the final state: 

U = K 00 -H 2 . (7.34) 

The conservation of U allows us to compute the kinetic part in the final state, 
K x , explicitly: 

K x = K - Hb + H 2 . (7.35) 

To obtain an identical asymptotic state for the HMF model, its kinetic energy 
must be equal to Koo and its magnetization to Ti. The energy for the HMF 
model is thus : 

Uhmf = K x + 1 (1 - Ti 2 ) . (7.36) 
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In addition to that energy, we impose that the height of the waterbag, /q, be 
the same for the two systems, as this quantity is an invariant of the dynamics. 

Explicitly, finding the set of AO and Ap amounts to solve the following 
equation: 

Ap 2 K sinA6» TT „„, 
-f + -(!-_) = %*, , (7.37) 

that transforms, making use of /o = 4Ae 1 Ap , into 

1 1,, sinA6> Tr 

96^ + 2 (1 -^ ) = W ■ (7 ' 38) 

We now present a comparison of the exact solution for the uncoupled pendula 
with a simulation of the HMF model. We consider a pendulum with an initial 
condition given by AO = 1.66 and Ap = 0.5. The resolution of the self-consistent 
problem for the exact solution gives b = 0.767. 

The initial condition for the HMF model, computed according to Eq.(7.38), 
is defined by AO = 0.523 and Ap = 1.390. The parameters for the simulation 
are N e = N p = 256 and At = 0.1. 
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Figure 7.16: Evolution of the magnetization M(t) in the HMF model in the 
course of time for the initial condition AO = 0.523, Ap = 1.390. The value 
predicted for the corresponding exact solution is displayed for reference. 

Figure 7.16 displays M(t) for the HMF simulation. M(t) begins by oscillat- 
ing strongly, then attains a value close to the one predicted by the theory. This 
similarity between the uncoupled pendula system and the HMF model comes a 
little bit as a surprise, given the difference between the two systems. We note 
that the value of M attained in the simulation is far from the original value 
Mq = 0.955, meaning that we can expect a better agreement when the final 
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and the initial magnetizations are closer, in a similar way to the results for 
Lynden-Bell's theory. 

For completeness, we also display the marginals for the two systems in 
Figs. 7.17 and 7.18, The general aspect of the curves is similar but a quan- 
titative agreement is not present. 




Figure 7.17: (left panel) Comparison Figure 7.18: (right panel) Same as 
of p{0) for a simulation of the HMF Fig. 7.18, but for <p(p). 
model and the exact solution for the 
corresponding system of pendula. See 
text for details. 



7.3.5 Occurrence of a phase transition 

Equation (7.32) always admit a solution where % = b = 0. It may possess a 
solution where H and b are finite. The passage from a region of parameters AO 
and Ap where two solutions exist to one where only the homogeneous solution 
is present is similar to a phase transition. This transition is found for instance 
in Ref. (Leoncini et al., 2009). 

We illustrate this behavior in Fig. 7.19 Above a threshold in Ap, H = b 
drops to 0. 



Ill 



0.40 




Figure 7.19: Representation of the solution of Eq. (7.32) for A8 = tt in function 
of Ap. We observe a transition from H = b finite to H = b = 0. 
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7.3.6 Discussion 



In this chapter, we have proposed a Vlasov equation for a system similar to the 
Hamiltonian Mean-Field (HMF) model but without interactions. We elaborated 
Lyndon-Bell's theory for that system and found a remarkable agreement for the 
predicted value of the bunching, in a range of parameters. 

With the help of the energy probability distribution functions (energy PDF) , 
we found that a factor explaining the effectiveness of the theory is the similarity 
between the initial energy PDF and the one predicted by the theory. This 
finding is an important addition to the understanding of violent relaxation in 
which incomplete relaxation is observed, as it confirms the adequacy of Lynden- 
Bell's theory for a wider range of applications than expected. 

Then we analyzed an exact solution of the asymptotic evolution of the Vlasov 
equation, under an crgodic-like hypothesis. This solution shows some similarity 
with the asymptotic state of the HMF model, especially with respect to the 
value of the bunching. Solving a self-consistency relation, we also displayed a 
transition between two regions of parameters in which the relation possesses 
cither one or two solutions. Understanding this transition with more detail is 
of great importance to the understanding of transitions in related mean-fields 
models such as the HMF model and the Free-Electron Laser. 
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Chapter 8 

Conclusions and 
perspectives 

We performed in this thesis the investigation of toy models in which the in- 
teraction potential is long ranged, namely the Hamiltonian Mean-Field model 
(HMF) and the Colson-Bonifacio model for the free-electron laser (FEL). Their 
dynamical description in the thermodynamic limit is cast with the help of the 
Vlasov equation. 

We reported the use of a numerical procedure to compute explicitly the 
evolution of these models. This represents the first detailed numerical com- 
putation for the Vlasov equation for simple models of long-range interactions. 
The present work puts emphasis on the interest of these computations for both 
domains of research : the numerical resolution of the Vlasov equation and the 
theoretical study of the Vlasov equation for models with long-range interactions. 

The occurrence of quasi-stationary states (QSS) and out-of-equilibrium phase 
transitions presents a conceptual challenge to the understanding of kinetic the- 
ory. At the present time, the statistical theory of Lynden-Bell provides a correct 
predictive and interpretative framework for the Vlasov dynamics. It has been 
reported to describe successfully the asymptotic dynamics of the HMF model 
and of the FEL in a number of situations. 

The study of the FEL via numerical simulations of the Vlasov equation 
provides a detailed description of the out-of-equilibrium dynamics. The intensity 
of the laser field produced by the FEL device takes significantly different values 
depending on the initial condition. Specializing on waterbag initial conditions, 
we described the out-of-equilibrium phase transition occurring between these 
regimes. We interpreted further this transition in term of Lynden-Bell's theory, 
but acknowledge that in the case of the FEL it does not yet provides a fully 
predictive theory. A direct observation of the structures in phase space allowed 
us to understand a shortcoming of the theory, as dynamical features cannot be 
understood within its framework. An original application of Vlasov simulations 
coupled to Poincare section of discrete particles is introduced and is expected 
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to yield accurate results on the phase space properties of the dynamics. The 
occurrence of pcriod-2 and of higher periods orbits is anticipated. 

We then turned to the investigation of the origin of relaxation in Vlasov dy- 
namics for the HMF model. The use of an analogy with fluid dynamics allows to 
quantify the deformations experienced by an initial waterbag : the evolution in 
time of the perimeter of the fluid gives information on the microscopic dynamics 
in phase space. The occurrence of stretching and folding structures, in relation 
with the filamentation, provides a new point of view on the collisionless evolu- 
tion and its approach to equilibrium. We performed the numerical computation 
of the perimeter in a situation and found it to be exponential. The high level 
of details required to compute the perimeter required the use of high accuracy 
simulations. 

The analogy between the HMF model and a set of uncoupled pendula is 
cast in the last chapter. We introduce the Vlasov equation corresponding to 
that system. Lynden-Bell's theory is developed and investigated in this context 
and is shown to predict an out-of-equilibrium phase transition under a self- 
consistency hypothesis. The results of Lynden-Bell's theory are encouraging, 
which is surprising for a system in which there is no interaction, hence no source 
of relaxation. A careful analysis reveals that the theory predicts good results 
when the initial condition presents an energy distribution function close to the 
one predicted by Lynden-Bell's theory. 

An exact computation of the asymptotic dynamics for this Vlasov equation 
is then performed. It matches exactly the simulation results up to a small 
filamentary structure. This system, under the self-consistency hypothesis, also 
displays a phase transition. A comparison is drawn with the HMF model, and 
we expect that the full understanding of the exact system will yield more insight 
on the out-of-equilibrium phase transition occurring in the HMF model or in 
the FEL. 

In the light of the results presented in this thesis, a number of questions re- 
main open. The exact characterization of QSSs is yet to be achieved. Especially, 
the dynamical features appearing in theses QSSs are not interpreted in terms 
of kinetic theory. The fine structures found in phase space have an impact on 
the macroscopic dynamics that needs to be investigated. 

The detailed presentation of numerical computations calls for comparisons 
with other methods with the hope that theoretical investigation on toy models 
may contribute to the progress of simulational methods. 

Finally, we would like to mention a number of encouraging results that point 
to possible experimental realizations of systems with long-range interactions. 

The out-of-equilibrium phase transition occurring in the Free-Electron Laser 
has the potential to be realized experimentally. A model for collective atomic 
recoil lasing (CARL) possesses a mathematical description close to the one of 
the FEL (Bonifacio and Salvo, 1994). Experiments on the CARL may offer more 
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flexibility and provide a new opportunity to investigate the effects of long-range 
interactions (Slama et al., 2008; Campa et al., 2009). 

The dipolar interaction in spin structures is also a candidate towards the 
experimental realization of a system with long-range interaction. The authors 
of Ref. (Campa et al., 2007) propose a setup in which the statistical mechan- 
ics of the system can display ensemble inequivalence, negative specific heat or 
temperature jumps. 
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Appendix A 

Stability of the 
homogeneous waterbag 

This Appendix is devoted to reviewing the stability condition of system (5.2) 
for an initial homogeneous distribution of the waterbag type. The derivation 
follows from a straightforward linear analysis which can be found for instance 
in (Elskens and Escande, 2003). We shall hereafter make reference to the cal- 
culation detailed in (Bachelard and Fanelli, 2010). Let us start by assum- 
ing a general equilibrium setting where the spatial distribution is homogeneous 
(A x = A y = 0) and / = fo(p), i.e. a generic function of the variable p. Then one 
can linearize around the equilibrium and eventually derive an explicit solution 
which holds for a relatively short time. To this end we write: 

f(6,p,t) = fo(p) + fi(6,p,t), 
A x (t) =X 1 (t) and A y (t)=Y 1 (t) . (A.l) 

where the quantities labeled with the index 1 stand for the linear perturbation. 
Introducing in system (5.2) and retaining the lowest order yields: 

(<9 2 + pd )h - M x i cos ~Yi sin 9) =0 (A.2) 

/TT Z'+OO 
d9 dp /i cos 9 =0 (A.3) 

-TT J — OO 

/TT r-\-oo 
d9 dp A sin0 +^ =0 (A.4) 
-TT J — OO 

where we introduced rj(p) = d p fo(p). The above linear system admits a solution 
in terms of normal modes: 

fi(6,p,z) = F^e^-^ +F*(p)e-« e -^^ (A.5) 
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Xi(z) = X 1 e~ iuz + XI e iu "~ z (A.6) 
Y x {z) = iYxe-™* - iY? e iu "~ z . (A.7) 

where the symbol * refers to the complex conjugate and in general u G C. 
Making use of the above ansatz in the linearized system of equations returns 
the following consistency equation: 

u = [ + °°dp^- (A.8) 

J-oo P-W 

often referred to as the dispersion relation. To determine whether a given dis- 
tribution fo (p) is stable or unstable, one can solve the above dispersion relation 
and estimate the sign of the imaginary part of uj. Depending on the sign the 
field grows exponentially (instability) or oscillates indefinitely (stability). If the 
selected initial condition is parametrized via an adjustable parameter, one can 
then calculate the corresponding theshold value which discriminates between 
stable and unstable regimes. This analytical procedure can be persecuted in 
simple cases, as the one addressed in this appendix (the waterbag). For more 
complicated situations one can resort to the celebrated Nyquist method, first 
introduced in plasma physics (Nyquist, 1932) (see also (Chavanis and Delfini, 
2009)). For the case at hand, fo{p) takes the form: 

fo(p) = ^^[&(p + Ap)-e(p-A P )} (A.9) 

where Q stands for the Heaviside function. Inserting (A.9) into (A.8), carry- 
ing out the integral explicitly and looking for the value of Ap which sets the 
transition between complex and real w, leads to Ap ~ 1.37 or equivalently 
e = (Apf/Q ~ 0.315. 
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Appendix B 

Routines 



We review in this appendix technical aspects of the computations presented in 
this thesis. 

B.l The simulation program 

The simulation program used to perform the Vlasov simulations has been writ- 
ten in the course of the thesis. It makes use of modern programming concept 
applied to the Fortran 90 language in the spirit of Ref. (Decyk et al., 1998). 

The program is based on a Fortran 90 module that contains the description 
of a grid in phase space. The main element of the module is the derived type 
"grid" that contains the basic physical values needed to perform the simulation. 
The derived type grid is shown in listing B.l . A number of subroutines makes use 
of "grid" and perform operations such as setting an initial condition, compute 
the marginal distributions, perform the advection steps in both directions of 
phase space on the basis of spline interpolation, . . . 

Specific modules for the Hamiltonian Mean-Field (HMF) model and the 
Free-Electron Laser combine "grid" with additional model-specific variables (the 
magnetization for the HMF model for instance). The Fortran 90 module con- 
tains in addition subroutines that write the data in a systematic way in hdf5 files. 
hdf5 is a structured file format available at http://www.hdfgroup.org/HDF5/ . 

B.2 Simulated annealing 

The simulated annealing algorithm (SA) is a general optimization algorithm 
(Kirkpatrick et al., 1983). It consists of a random sampling in the space of 
parameters in which each moves is accepted or rejected on the basis of a cost 
function. It is used in this thesis to find "candidate roots" to the set of equations 
given by Lynden-Bell's theory. 

The SA algorithm accepts any move that decreases the cost function. For 
the moves that increase the cost function, the acception is based on the following 



121 



Listing B.l: Derived data type used in the Fortran 90 simulation program, 
type grid 

integer :: Nx, Nv / The number of grid points 

in the space and the 

velocity directions 
double precision :: xmin , xmax, vmin , vmax 
/ / The boundaries of phase 

space 

double precision : : dx , dv / The step of the grid in 

the space and the velocity 
directions 

double precision : : DT / The time step 

double precision, allocatable :: f(: ,:) 

/ / The array containing f. 

double precision , allocatable : : f 2 ( : , : ) , g ( : , : ) 

/ / The arrays containing the 

second derivatives of f 
needed to perform the 
spline interpolation (f2) 
and a copy of f (g) . 
double precision, allocatable :: rho ( : ) , phi(:) 
/ / The arrays containing the 

marginals in the x and in 
v . 

end type grid 



criterion : 

— Acost_function 
exp >£ (B.l) 

where cost_function is the cost function, T is the effective temperature and 
< £ < 1 is a random number. The structure of such a program is given in 
listing B.2. 

Performing a series of SA runs allows one to find candidates for the solution 
of the set of equations defined by Lynden-Bell's theory. An example is given in 
Fig. B.l : The minimum of 16 runs for the same parameters is displayed. This 
minimum decreases in the course of the iterations. The programmed tempera- 
ture is defined as T(i) = exp 2 - 5 . 1 . 
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Listing B.2: Principle of the simulated annealing algorithm, 
do i=l,n_steps 



T = T_program(i) 

oldx = x 

call move(x, T) 

cost = cost_function(x) 



/ T ^program is the 
programmed temperature . 
! We keep the old parameter 

in case of rejection . 
! We perform a random move 
of the parameter x. 
! co st_f unction is a cost 
function defined by the 
user 

if ( cost . leq . oldcost ) then 

.' We accept the move if the cost decreases . 

oldcost = cost 
else 

if (exp( — (cost — oldcost) / T) . gt . randomjiumber ( ) ) 
then 

oldcost = cost 
else 

/ else we reject the move 
x = oldx 
end i f 
end i f 
end do 




2000 4000 6000 8000 10000 

i 



Figure B.l: Evolution of the minimum cost in a series of 16 runs of simulated 
annealing. The computation is performed for the set of equations (7.14) with 
parameters A8 — 2.5 and Ap = 0.4. 



123 



124 



Bibliography 



Antoni, M. and Ruffo, S. (1995). Clustering and relaxation in hamiltonian 
long-range dynamics. Phys. Rev. E, 52(3):2361. 

Antoniazzi, A., Califano, F., Fanelli, D., and Ruffo, S. (2007a). Exploring the 
thermodynamic limit of hamiltonian models: Convergence to the vlasov 
equation. Phys. Rev. Lett., 98(15):150602. 

Antoniazzi, A., De Ninno, G., Fanelli, D., Guarino, A., and Ruffo, S. (2005). 
Wave-particle interaction: from plasma physics to the free-electron laser. 
Journal of Physics : Conference Series, 7:143. 

Antoniazzi, A., Fanelli, D., Barre, J., Chavanis, P.-H., Dauxois, T., and Ruffo, 
S. (2007b). Maximum entropy principle explains quasistationary states 
in systems with long-range interactions: The example of the hamiltonian 
mean-field model. Phys. Rev. E, 75(1):011112. 

Antoniazzi, A., Fanelli, D., Ruffo, S., and Yamaguchi, Y. Y. (2007c). Nonequi- 
librium tricritical point in a system with long-range interactions. Phys. 
Rev. Lett, 99(4):040601. 

Arbcr, T. and Vann, R. G. L. (2002). A critical comparison of eulerian-grid- 
based vlasov solvers. J. Corny. Phys., 180(1):339. 

Bachelard, R., Chandrc, C, Fanelli, D., Leoncini, X., and Ruffo, S. (2008). 
Abundance of regular orbits and nonequilibrium phase transitions in 
the thermodynamic limit for long-range systems. Phys. Rev. Lett., 
101(26):260603. 

Bachelard, R. and Fanelli, D. (2010). Short-time dynamics in the presence of 
wave-particles interactions: A perturbative approach. Commun. Nonlin- 
ear. Sci. Numer. Simulat., 15:40. 

Balcscu, R. (1997). Statistical Dynamics - Matter out of Equilibrium. Imperial 
College Press, Imperial College - London. 

Barre, J. (2003). Mecanique Statistique Et Dynamique Hors Equilibre de 
Systemes avec Interactions a Longue Portee. PhD thesis, Ecole Normale 
Superieure de Lyon. 



125 



Barrc, J., Bouchet, F., Dauxois, T., and Ruffo, S. (2005). Large deviation 
techniques applied to systems with long-range interactions. J. Stat. Phys., 
119(3-4):677. 

Barrc, J., Bouchet, F., Dauxois, T., Ruffo, S., and Yamaguchi, Y. Y. (2006). 
The vlasov equation and the hamiltonian mean- field model. Physica A, 
365:177. 

Barre, J., Dauxois, T., De Ninno, G., Fanelli, D., and Ruffo, S. (2004). Sta- 
tistical theory of high-gain free-electron laser saturation. Phys. Rev. E, 
69(4):045501. 

Barrc, J., Mukamcl, D., and Ruffo, S. (2001). Inequivalence of ensembles in a 
system with long-range interactions. Phys. Rev. Lett., 87(3):030601. 

Bcnedctti, C, Rambaldi, S., and Turchetti, G. (2006). Relaxation to boltzmann 
equilibrium of 2d coulomb oscillators. Physica A, 364:197. 

Bonifacio, R., Casagrande, F., Cerchioni, G., Souza, L. D. S., Pierini, P., and 
Piovella, N. (1990). Physics of the high-gain fel and superradiance. Rivista 
del Nuovo Cimento, 13(9):1. 

Bonifacio, R., Casagrande, F., and Pellegrini, C. (1987). Hamiltonian model of 
a free electron laser. Opt. Commun., 61(1):55. 

Bonifacio, R., Casagrande, F., and Souza, L. D. S. (1986). Collective variable 
description of a free-electron laser. Phys. Rev. A, 33(4):2836(R). 

Bonifacio, R., Pellegrini, C, and Narducci, L. M. (1984). Collective instabilities 
and high-gain regime in a free electron laser. Opt. Commun., 50(6):373. 

Bonifacio, R. and Salvo, L. D. (1994). Exponential gain and self-bunching in a 
collective atomic recoil laser. Phys. Rev. A, 50(2):1716. 

Bourke, P. (1987). Conrec - a contouring subroutine. Byte Magazine. 

Braun, W. and Hepp, K. (1977). The vlasov dynamics and its fluctuations in the 
1/ n limit of interacting classical particles. Commun. Math. Phys., 56:101. 
(c) 1977: Springer- Verlag. 

Brinkmann, R. (2006). The european xfel project. In Proceedings of the FEL 
conference 2006, Berlin, Germany. 

Califano, F., Galeotti, L., and Mangeney, A. (2006). The vlasov-poisson model 
and the validity of a numerical approach. Phys. Plasmas, 13(8):082102. 

Campa, A., Dauxois, T., and Ruffo, S. (2009). Statistical mechanics and dy- 
namics of solvable models with long-range interactions. Phys. Rep., 480:57. 



126 



Campa, A., Giansanti, A., Morigi, G., and Sylos Labini, F., editors (2008). 
Dynamics and Thermodynamics of Systems with Long Range Interactions: 
Theory and Experiments, Assist, Italy J^-8 July 2007, volume 970 of AIP 
Conference Proceedings. American Institute of Physics. 

Campa, A., Khomeriki, R., Mukamel, D., and Ruffo, S. (2007). Long-range 
effects in layered spin structures. Phys. Rev. B, 76:64415. 

Canosa, J., Gazdag, J., and Fromm, J. (1974). The recurrence of the initial 
state in the numerical solution of the vlasov equation. J. Comp. Phys., 
15(1):34. 

Chandrasckhar, S. (1942). Principles of stellar dynamics. The University of 
Chicago Press, Chicago, Illinois. 

Chavanis, P.-H. (2006a). Hamiltonian and brownian systems with long-range in- 
teractions: I. statistical equilibrium states and correlation functions. Phys- 
ica A, 361:55. 

Chavanis, P.-H. (2006b). Hamiltonian and brownian systems with long-range 
interactions: II. kinetic equations and stability analysis. Physica A, 361:81. 

Chavanis, P.-H. (2006c). Lynden-bell and tsallis distributions for the hmf model. 
Eur. Phys. J. B, 53(4):487. 

Chavanis, P.-H. (2006d). Quasi-stationary states and incomplete violent relax- 
ation in systems with long-range interactions. Physica A, 365:102. 

Chavanis, P.-H. (2008a). Hamiltonian and brownian systems with long-range 
interactions: III. the bbgky hierarchy for spatially inhomogeneous systems. 
Physica A, 387:787. 

Chavanis, P.-H. (2008b). Hamiltonian and brownian systems with long-range in- 
teractions: IV. general kinetic equations from the quasilinear theory. Phys- 
ica A, 387:1504. Elsevier Ltd. 

Chavanis, P.-H. (2008c). Hamiltonian and brownian systems with long-range 
interactions: V. stochastic kinetic equations and theory of fluctuations. 
Physica A, 387:5716. 

Chavanis, P.-H. and Dclfini, L. (2009). Dynamical stability of systems with long- 
range interactions: application of the nyquist method to the hmf model. 
Eur. Phys. J. B, 69(3):389. 

Chavanis, P.-H., Sommeria, J., and Robert, R. (1996). Statistical mechanics 
of two-dimensional vortices and collisionless stellar systems. Astrophys. J., 
471(l):385-399. 

Cheng, C. Z. and Knorr, G. (1976). The integration of the vlasov equation in 
configuration space. J. Comp. Phys., 22:330. 



127 



Chomaz, P. and Gulminelli, F. (2002). Phase transitions in finite systems. In 
Dauxois, T., Ruffo, S., Arimondo, E., and Wilkens, M., editors, Dynamics 
and Thermodynamics of Systems With Long Range Interactions, volume 
602 of Lecture Notes in Physics, page 68. Springer- Ver lag. 

Colson, W. B. (1976). Theory of a free electron laser. Physics Letters A, 59:187. 

Crouseilles, N., Latu, G., and Sonnendrucker, E. (2007). Hermite spline inter- 
polation on patches for parallelly solving the vlasov-poisson equation. Int. 
J. Appl. Math. Comput. ScL, 17(3):335. 

Curbis, F., Antoniazzi, A., De Ninno, G., and Fanelli, D. (2007). Maximum 
entropy principle and coherent harmonic generation using a single-pass free- 
electron laser. Eur. Phys. J. B, 59:527. 

Dauxois, T., Ruffo, S., Arimondo, E., and Wilkens, M., editors (2002). Dy- 
namics and Thermodynamics of Systems With Long Range Interactions, 
volume 602 of Lecture Notes in Physics. Springer- Verlag. 

Dauxois, T., Ruffo, S., and Cugliandolo, L. (2008). Les houches summer school 
on long-range interacting systems. 

Dawson, J. M. (1983). Particle simulation of plasmas. Rev. Mod. Phys., 55:403. 

de Buyl, P. (2009). Numerical resolution of the vlasov equation for the hamil- 
tonian mean-field model. Commun. Nonlinear. Sci. Numer. Simulat. 

de Buyl, P., Fanelli, D., Bachclard, R., and De Ninno, G. (2009). Out-of- 
equilibrium mean-field dynamics of a model for wave-particle interaction. 
Phys. Rev. ST Accel. Beams, 12(6):060704. 

de Buyl, P., Mukamcl, D., and Ruffo, S. (2005). Ensemble inequivalence in a 
xy model with long-range interactions. AIP Conf. Proc, 800:533. 

De Ninno, G., Allaria, E., Coreno, M., Chowdhury, S., Curbis, F., Danailov, 
M. B., Diviacco, B., Ferianis, M., Karantzoulis, E., Longhi, E. C, Pinayev, 
I. V., Spezzani, C, Trovo, M., and Litvinenko, V. N. (2008). Self-induced 
harmonic generation in a storage-ring free-electron laser. Phys. Rev. Lett., 
100(10):104801. 

Decyk, V. K., Norton, C. D., and Szymanski, B. K. (1998). How to support 
inheritance and run-time polymorphism in fortran 90. Comput. Phys. Com- 
mun., 115:9. 

del Castillo-Negrete, D. (2002). Dynamics and self-consistent chaos in a mean- 
field hamiltonian model. In Dauxois, T., Ruffo, S., Arimondo, E., and 
Wilkens, M., editors, Dynamics and Thermodynamics of Systems With 
Long Range Interactions, volume 602 of Lecture Notes in Physics, page 
407. Springer- Verlag. 



128 



Doyuran, A., Babzien, M., Shaftan, T., Yu, L. H., Dimauro, L. F., Ben-Zvi, L, 
Biedron, S. G., Graves, W., Johnson, E., Krinsky, S., Malone, R., Pogorel- 
sky, I., Skaritka, J., Rakowsky, G., Wang, X. J., Woodle, M., Yakimenko, 
V., Jagger, J., Sajaev, V., and Vasserman, I. (2001). Characterization of a 
high-gain harmonic-generation free-electron laser at saturation. Phys. Rev. 
Lett, 86(26):5902. 

Drewsen, M., Mortensen, A., Nielsen, E., and Matthey, T. (2008). Strongly 
correlated ion coulomb systems. In Campa, A., Giansanti, A., Morigi, G., 
and Labini, F. S., editors, Dynamics and Thermodynamics of Systems with 
Long-Range Interactions: Theory and Experiments, volume 970, page 295. 
AIP. 

Driscoll, C. F. and Fine, K. S. (1990). Experiments on vortex dynamics in pure 
electron plasmas. Phys. Fluids B, 2:1359. 

Elskcns, Y. and Escande, D. (2003). Microscopic Dynamics of Plasmas and 
Chaos. Institute of Physics Publishing. 

Filbet, F., Sonnendrucker, E., and Bcrtrand, P. (2001). Conservative numerical 
schemes for the vlasov equation. J. Comp. Phys., 172:166. 

Firpo, M.-C. (2009). Unveiling the nature of out-of-equilibrium phase transitions 
in a system with long-range interactions. Europhys. Lett., 88:30010. 

Funato, Y., Makino, J., and Ebisuzaki, T. (1992). Violent relaxation is not a 
relaxation process. Publ. Astron. Soc. Japan, 44:613. 

Galeotti, L. and Califano, F. (2005). Asymptotic evolution of weakly collisional 
vlasov-poisson plasmas. Phys. Rev. Lett., 95:15002. 

Gutnic, M., Haefele, M., Paun, I., and Sonnendrucker, E. (2004). Vlasov simula- 
tions on an adaptive phase-space grid. Comput. Phys. Commun., 164:214. 

Hockney, R. W. and Eastwood, J. W. (1988). Computer Simulations Using 
Particles. Adam Hilger, Bristol and New- York. 

Hohl, F. (1968). Theory and results on collective and collisional effects for a 
one-dimensional self-gravitating system. NASA Technical Report, TR:R- 
289. 

Inagaki, S. and Konishi, T. (1993). Dynamical stability of a simple model similar 
to self-gravitating systems. Publ. Astron. Soc. Japan, 45:733. 

Jain, K., Bouchet, F., and Mukamel, D. (2007). Relaxation times of unstable 
states in systems with long range interactions. J. Stat. Mech., 11:11008. 

Kirkpatrick, S., Gclatt, C. D., and Vecchi, M. P. (1983). Optimization by 
simulated annealing. Science, 220:671. 



129 



Latora, V., Rapisarda, A., and Ruffo, S. (1998). Lyapunov instability and finite 
size effects in a system with long-range forces. Phys. Rev. Lett., 80(4):692. 

Leoncini, X., Berg, T. V. D., and Fanelli, D. (2009). Out of equilibrium solutions 
in the xy-hamiltonian mean- field model. Europhys. Lett., 86:20002. 

Levin, Y., Pakter, R., and Rizzato, F. B. (2008a). Collisionless relaxation in 
gravitational systems: From violent relaxation to gravothermal collapse. 
Phys. Rev. E, 78(2):021130. 

Levin, Y., Pakter, R., and Teles, T. N. (2008b). Collisionless relaxation in 
non-neutral plasmas. Phys. Rev. Lett., 100(4):40604. 

Luwel, M. and Severne, G. (1985). Collisionless mixing in 1-dimensional grav- 
itational systems initially in a stationary waterbag configuration. A&A, 
152:305. 

Luwel, M., Severne, C, and Rousseeuw, P. J. (1984). Numerical study of the 
relaxation of one-dimensional gravitational systems. Astrophys. Space ScL, 
100:261. 

Lynden-Bell, D. (1967). Statistical mechanics of violent relaxation in stellar 
systems. Mon. Not. R. Astron. Soc. 

Menotti, C, Lewenstein, M., Lahaye, T., and Pfau, T. (2008). Dipolar interac- 
tion in ultra-cold atomic gases. In Campa, A., Giansanti, A., Morigi, G., 
and Labini, F. S., editors, Dynamics and Thermodynamics of Systems with 
Long-Range Interactions: Theory and Experiments, volume 970, page 332. 
AIP. 

Mineau, P., Feix, M. R., and Rouet, J. L. (1990). Numerical simulations of vio- 
lent relaxation and formation of phase space holes in gravitational systems. 
A&A, page 6. 

Morita, H. and Kaneko, K. (2006). Collective oscillation in a hamiltonian sys- 
tem. Phys. Rev. Lett., 96(5):050602. 

Nyquist, H. (1932). Regeneration theory. Bell Systems Technical Journal, 
11:126. 

Ottino, J. M. (1989). The kinematics of mixing: stretching, chaos and transport. 
Cambridge University Press, Cambridge, England. 

Pichon, C. (1994). Dynamics of self-gravitating disks. PhD thesis, University of 
Cambridge. 

Press, W. H., Tcukolsky, S. A., Vetterling, W. T., and Flannery, B. P. (1996). 
Numerical Recipes in Fortran 90. Cambridge University Press, 2nd edition. 



130 



Ruffo, S. (1994). Hamiltonian dynamics and phase transitions. In Benkadda, S., 
Elskens, Y., and Doveil, F., editors, Transport and Plasma Physics. World 
Scientific, Singapore. 

Severne, G. and Luwel, M. (1986). Violent relaxation and mixing in non-uniform 
one-dimensional gravitational systems. Astrophys. Space Sci., 122:299. 

Sguanci, L., Gross, D., and Ruffo, S. (2005). Apparent fractal dimensions in the 
HMF model. Transport Theory and Statistical Physics, 34(3-5) :431. 

Shoucri, M. (2008). Eulerian codes for the numerical solution of the vlasov 
equation. Commun. Nonlinear. Sci. Numer. Simulat., 13(1):174. 

Slama, S., Krenz, G., Bux, S., Zimmermann, C, and Courteille, P. W. (2008). 
Collective atomic recoil lasing and superradiant rayleigh scattering in a 
high-Q ring cavity. In Campa, A., Giansanti, A., Morigi, G., and Labini, 
F. S., editors, Dynamics and Thermodynamics of Systems with Long-Range 
Interactions: Theory and Experiments, volume 970, page 319. AIP. 

Sonncndrucker, E., Roche, J., Bertrand, P., and Ghizzo, A. (1999). The semi- 
lagrangian method for the numerical resolution of the vlasov equation. J. 
Comp. Phys., 149:201. 

Staniforth, A. and Cote, J. (1991). Semi-lagrangian integration schemes for 
atmospheric models - a review. Mon. Wea. Rev., 119:2206. 

Watanabe, T.-H. and Sugama, H. (2005). Vlasov and drift kinetic simulation 
methods based on the symplectic integrator. Transport Theor. Stat. Phys., 
34:287. 

Yamaguchi, Y. Y. (2008). One-dimensional self-gravitating sheet model and 
lynden-bell statistics. Phys. Rev. E, 78(4):041114. 

Yamaguchi, Y. Y., Barre, J., Bouchet, F., Dauxois, T., and Ruffo, S. (2004). 
Stability criteria of the vlasov equation and quasi-stationary states of the 
hmf model. Physica A, 337(l-2):36. 

Yamashiro, T., Gouda, N., and Sakagami, M. (1992). Origin of core-halo 
structure in one-dimensional self-gravitating system. Prog. Theor. Phys., 
88(2):269. 

Zwillinger, D. (1998). Handbook of Differential Equations. Academic Press, San 
Diego - CA, 3rd edition. 



131 



